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PREFACE 

This Memorandum describes a mathematical simulation model of sea 
ice growing and dissipating with the seasons and interacting with the 
atmosphere and the ocean. The model was developed by G. A. Maykut, 
performing research at the University of Washington under supervision 
of Rand consultant N. Untersteiner. It was a product of an arctic 
research program sponsored by the Office of Naval Research, which be- 
cause of its immediate application to Rand weather studies has con- 
curred with the publication of the work by Rand. 

The subject is interrelated with a continuing line of Rand inves- 
tigations into climatic change, its nature, its causes, and the possi- 
bilities for influencing climate (whether deliberately or inadvertent- 
ly) on a large scale. The general aim of the program is to understand 
the dynamics of the thermally forced motion of the global atmosphere 
and ocean and to simulate this behavior using mathematical models and 
high-speed computers. Related work by other Rand staff members and 
consultants includes the development of simulation models of the atmos- 
phere (Y. Mintz and A. Arakawa) and of the wind-driven ocean (W. L. 
Gates), studies of large-scale variations in thermal forcing of the 
global system (J. Bjerknes, J. 0. Fletcher, R. E. Huschke, S. Oleni- 
coff, and R. R. Rapp) , and the development of simulation models of con- 
vective processes which heat the atmosphere and produce rain (L. R. 
Koenig and F. W. Murray). 

Variations in ice extent on the sea are associated with variations 
in thermal forcing of the global system. Previous Rand reports explor- 
ing this interaction include: The Heat Budget of the Arctic Basin and 
its Relation to Climate . R-444, 1965, Proceedings of the Symposium on 
the Arctic Heat Budget and Atmospheric Circulation , RM-5233, 1966, and 
Ice Extent on the Southern Ocean and its Relation to World Climate , 
RM-5793, 1969. This Memorandum represents one more step toward under- 
standing and simulating the dynamics of the global system. 



SUMMARY 



A one-dimensional thermodynamic model of sea ice is presented 
which includes the effects of snow cover, ice salinity, and internal 
heating due to penetration of solar radiation. Surface-energy balances 
determine rates of ablation and accretion; diffusion equations govern 
heat transport within the ice and snow. The incoming radiative and 
turbulent fluxes, oceanic heat flux, ice salinity, snow accumulation, 
and surface albedo are specified as functions of time. Starting from 
an arbitrary initial condition, the model is integrated numerically 
until annual equilibrium patterns of temperature and thickness are 
achieved. The model neglects the effects of mechanical stresses on 
the ice and the influence of surface conditions on the turbulent fluxes. 
Leads and melt ponds are considered only insofar as they indirectly af- 
fect the energy fluxes over the bare ice. The model is applied to the 
Central Arctic. 

The input data for the initial test of the model were drawn 
largely from J. 0. Fletcher, The Heat Budget of the Arctic Basin and 
Its Relation to Climate , The Rand Corp., Santa Monica, (1965). Values 
predicted by the model for the average ice thickness (288 cm) , amount 
of surface ablation (40 cm) , and the temperature field all agree closely 
with field observations. To define the response of sea ice to changes 
of its environment, various sets of values for the input parameters were 
specified, and additional solutions obtained. The roles of oceanic heat 
flux, snow cover, internal heating, albedo, incoming energy fluxes, and 
ice salinity were examined. Under present conditions, the ocean must 
supply 1 to 2 kcal/cm 2 year to the ice; an additional 4 kcal/cm 2 year 
would cause the ice to vanish. Analysis of these results has also shown 
( that the response time of sea ice for a given set of parameters is de- 
termined primarily by the final equilibrium thickness, rather than the 
initial thickness. Annual snow depths less than 70 cm have little ef- 
fect on equilibrium thickness because of mutually compensating changes 
in ice temperature, ice ablation, and internal heating; snow depths 
greater than 70 cm result in much thicker ice. Comparison of observed 
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and calculated temperature profiles indicate that 2.0 to 2.5 kcal/cm 2 year 
of the short-wave radiation penetrates the ice and contributes to inter- 
nal heating, but as a larger fraction is allowed to penetrate the ice, 
less energy is available for melting at the surface and the average 
thickness increases. An increase in the amount of either incoming long- 
wave or short-wave radiation by 10 percent results in a 50 percent de- 
crease in ice thickness. Although the turbulent fluxes are small in com- 
parison to the radiative fluxes, they are significant in terms of the 
net energy balance. If their gum is reduced from -0.5 kcal/cm 2 year to 
zero, a 70 percent decrease in thickness occurs. Average ice albedos 
under 0.50 cause the ice to vanish in a few years. 

The present model indicates that, of the various techniques pro- 
posed for large-scale ice removal, those involving modification of the 
surface would be the most feasible and effective. It is suggested that, 
in future field projects, priority be given to studies investigating the 
interaction of leads and melt ponds with the ice. 
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I . INTRODUCTION 
1.1 SEA ICE AND THE POLAR REGIONS 

1 The polar regions, as the major heat sinks of the Earth, influence 
the general circulation patterns in both hemispheres. The heating at 
high and low latitudes is unequal, and the consequent general horizon- 
tal temperature gradient from the polar regions to the lower latitudes 
forces heat into the polar regions. In the Northern Hemisphere the 
annual heat transport into the region bounded by the Arctic Circle is 
approximately 75 kcal/cm^, averaged over the entire area. (This is 
roughly equivalent to the yearly total of incoming short-wave radia- 
tion reaching the surface of the Arctic.) Most of this advected heat 
is lost in space in the form of long-wave radiation, thereby helping 
to maintain the radiative equilibrium of the earth. 

•Heat also is carried into the Arctic Basin by ocean currents, 
but only about one-tenth as much as is carried by the atmosphere because 
of the relative slowness of oceanic transport and the relative inef- 
ficiency of the transfer of heat between ocean and atmosphere. Still, 
the oceanic heat flux is sufficient to affect the sea ice significantly, 
both in thickness and extent. 

Sea ice is a dominant feature of both polar oceans . The Arctic 

17 2 

Ocean (as defined by Zubov , 1943), covering some 10 cm , is about 
80 percent ice-covered during the winter and 60 percent in the late 
summer. The pack ice surrounding the Antarctic Continent has a maxi- 
mum area over one and a half times as great as in the Arctic. The 
annual variation in areal extent is about 85 percent ( Treshnikov , 1967) , 
contrasting with 25 percent in the Arctic. These broad changes make 
sea ice one of the most variable physical features on the surface of 
the earth. 

Because of differences in the configuration of oceans and land 
masses, sea ice plays different roles in the Arctic and the Antarctic. 
In the Arctic, with its ice-covered ocean, the extent of sea ice is 
critical to the energy balance over the entire region. Even small 
anomalies in the energy fluxes at the surface can greatly vary the 
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percentage of ocean that is ice-covered. In fact, a mere 22 kcal/cm 2 
surplus (only a third of the yearly energy flux into the Arctic) would 
obliterate the entire ice pack. The case of the Antarctic is strik- 
ingly different. There, the sea ice surrounding the continent, while 
greatly influencing the regional heat budget and the seasonal varia- 
tions in atmospheric heat advection, still probably affects the gen- 
eral features of atmospheric circulation only about as much as do the 
stable continental ice masses. The continental region, under its 2000 
to 4000 meters of ice is comparable in extent to the area covered by 
pack ice in the Arctic. Although antarctic sea ice could vanish with 
much the same energy input as that required for the Arctic, the masses 
of continental ice could be significantly affected by nothing less 
than a drastic climate change. Rapid removal of the southern pack ice 
would produce changes of unknown magnitude, but changes certainly ame- 
liorated by the stabilizing effect of the relatively permanent ice 
sheets, and no radical change in the general circulation would neces- 
sarily follow. 

For studying processes in the Northern Hemisphere, the Arctic 
appears to be of considerable concern, and the situation also is less 
complex, in the absence of large-scale continental ice. Considering 
further that the great annual variation in the extent of antarctic sea 
ice has made research difficult, the Arctic lends itself more readily 
to studies of sea ice and its large-scale effects. The following para- 
graphs will examine in more detail the role of sea ice in the Arctic. 

1.2 CLIMATIC CONSEQUENCES OF 

ICE COVERING THE ARCTIC SEA 

Sea ice modifies its environment in a number of ways. During the 
dark winter months, intense atmospheric cooling occurs; surface tem- 
peratures in the Arctic usually reach -30 to -40°C. The resulting in- 
crease in latitudinal temperature gradients produces an increase in 
the general circulation which then causes large quantities of warm air 
to be drawn into the Arctic from lower latitudes. Moreover, an ice 
cover transfers little heat from the ocean to the atmosphere, since 
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molecular conduction is the controlling mechanism. During the summer 
the pack ice is even more influential. Because of a relatively high 
albedo, only 30 to 40 percent of the incoming short-wave radiation is 
absorbed by an ice surface, as compared to 90 percent by open water. 
The large amount of radiation reflected contributes somewhat to atmos- 
pheric warming and a reduction in the intensity of atmospheric circu- 
lation. Heat storage in the Arctic Ocean is severly limited, while 
wind stress and wind mixing in the ocean is suppressed. Because little 
of the ocean surface is open for evaporation, arctic winters are char- 
acterized by light snowfall. Even during cloudy summers, precipitation 
remains low, being on the order of 10 to 15 g/cm^ year. 

The climatic consequences of ice covering the Arctic Ocean can 
best be understood through an examination of some possible results of 
an Arctic Ocean freed of ice. Fletcher (1965) , in an extensive eval- 
uation of existing literature dealing with the heat budget of the Arc- 
tic, visualizes an ice-free Arctic in roughly the following manner: 
in the summer the ocean would absorb and store up to 90 percent of all 
incoming solar radiation. This would imply a slightly cooler atmos- 
phere and hence a small increase in the intensity of atmospheric cir- 
culation. Greater amounts of evaporation would produce general cloud- 
iness over the area. During the winter the ocean would slowly release 
heat stored from the summer. Surface temperatures would be slightly 
above freezing in contrast to -35°C under present conditions. Advec- 
tion of heat from lower latitudes would be only slightly less than at 
present. An ice-free Arctic would thus be characterized by a constant, 
vigorous, year-round zonal flow with cool moist summers and warm moist 
winters. This situation should, according to Fletcher, induce heavy 
snowfall along coastlines with onshore flow and might possibly lead to 
accumulation of continental ice sheets. Quantitative results from a 
model of zonal temperature distribution (Rakipova , 1966) generally 
tend to corroborate Fletcher, although the model predicts a slight 
decrease in the intensity of both zonal and meridional circulation in 
the summer, as well as winter. In either case, the climatic changes 
in the more populous regions of the Northern Hemisphere could be dras- 
tic indeed. The primary difficulty in evaluating the magnitude of 
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any of these effects is that the complex interrelationship between 
the atmosphere and the oceans cannot yet be adequately described. 

1.3 RESPONSE OF ARCTIC SEA ICE TO ITS ENVIRONMENT 

The extent and the thickness distribution of the pack ice depend 
on conditions in the atmosphere and ocean. Laikhtman (1959) investi- 
gated the ways in which the atmosphere may affect the speed of surface 
ablation and concluded that air temperature and cloudiness are the two 
major factors controlling ice melt. Anomalies of cloudiness can cause 
great departures from the normal radiation balance, whereas higher air 
temperatures raise the sensible heat flux. He also suggests that wind 
speed would have little influence on summer melting because a larger 
downward sensible heat flux would be offset by a correspondingly larger 
evaporative heat loss. Budyko (1966) contends that the ice is so sen- 
sitive to temperature that a positive summer anomaly of only 4*C would 
cause the entire ice pack to melt in four years; a 2-degree anomaly 
would produce the same result in a few decades. However, the validity 
of these conclusions is open to some doubt. 

Although the winds exert little direct influence on surface abla- 
tion, they provide the principal mechanism through which ice movements 
occur. Figure 1 shows the average pressure patterns at high latitudes, 
superimposed over the mean ice drift. It may be seen that the ice mo- 
tions tend to parallel the isobars. Two important features arise from 
the motions of wind-driven ice: pressure ridges and open leads. Over- 
thrusting by ice flows can create zones where the ice thickness is up 
to an order of magnitude greater than that of the surrounding ice. 
Wittmann and Schule (1966) have found that 13 to 18 percent of the 
ice-covered portion of the Canadian Basin consists of pressure ridges, 
jlf pressure ridging proves to be this extensive throughout the Arctic 
Basin, existing estimates of the total amount of ice in the Arctic must 
be reviewed. Leads, especially during the winter, are somes of intense 
heat loss. Badgley (1961) found experimentally that the heat loss is 
at least two orders of magnitude larger from open leads than from per- 
ennial ice. Areas of open water thus stimulate rapid ice production. 




Fig. 1 - 



•- Annual mean patterns of surface pressure and ice movement 
in the Arctic ( Fletcher , 1965; data from Dunbar and 
Wittman , 1963). 
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The oceans influence arctic sea ice in only two principal ways: 
through turbulent heat transfer from the water to the ice and ice ex- 
port by winds and surface currents. The oceanic heat flux ultimately 
results from the advection of warmer water into the Arctic Ocean; the 
rate of heat transfer depends upon the magnitude of water exchange, 
the ocean temperatures, and the density stratification in the ocean. 
The primary exchange mechanism is the influx of warm Atlantic water 
through the Faroe/Shetland Channel and the outflux of cold arctic 
water via the East Greenland Current. A small amount of water enters 
by the Bering Strait and by continental runoff from rivers, but the 
heat transported by these sources is an order of magnitude less than 
that from the Atlantic. Nearly all the ice leaving the Arctic drifts 
through the Greenland/Spitsbergen Strait. The exact total of this ice 
export is in some doubt; many authors prefer the estimate of Zubov 
(1950) of 3000 km 3 per year. This figure represents roughly 10 per- 
cent of the total ice contained in the Arctic Basin. If the ice ex- 
ported each year were evenly distributed over the entire Polar Ocean, 
it would form a layer some 30 cm in thickness. Assuming that the pack 
ice is in near-equilibrium, average yearly ice production must lie 
near this value. The heat added to the region by the export of un- 
melted ice thus is about 2 kcal/cm 2 year. 

Because the water transfer is restricted and the ice cover is 
widespread, the structure of the Arctic Ocean is quite different from 
that of other oceans. As the ice insulates the ocean from the atmos- 
phere and reduces mechanical mixing, it is the primary determinant of 
vertical structure. The structure of the Arctic Ocean is horizontally 
quite uniform and characterized by an extremely stable density gradi- 
ent. However, during periods of rapid ice accretion and consequent 
brine expulsion by the ice, a shallow convective layer is formed. 
Vertical overturning of this unstable layer brings up warmer water 
from the lower layers, generating an increased heat flux toward the 
bottom of the ice. Therefore, analogously to the atmosphere, the 
ocean not only influences the extent and thickness of the ice, but 
also is itself modified by the ice. 
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1.4 SURVEY OF PREVIOUS RESEARCH 

The scientific literature dealing with the Arctic has multiplied 
many-fold since the pioneering voyages of the Fram (1893-1896) and the 
Maud (1918-1925). In 1937 the Soviets began a systematic study of the 
arctic interior by setting up the drifting station "North Pole 1 (NP-1)" 
on the pack ice. With the establishment of "Fletcher's Ice Island 
(T-3)" in 1952, the United States initiated its own long-range scien- 
tific program in the Arctic. Both the American and Soviet research 
efrorts continue. Recently data from the ice stations have been sup- 
plemented by observations from air-borne and satellite-borne instru- 
ments and submarine transits. The resulting mass of observational and 
experimental data has been summarized in several comprehensive surveys 
which are discussed in this section. 

By 1940 a coherent, although somewhat crude, picture of the Arctic 
and its physical environment had begun to emerge. N. N. Zubov in 1943 
published "Arctic Ice," a classic reference devoted entirely to arctic 
ice and oceanography. Zubov' s presentation of the geophysical aspects 
of the Arctic was, at the time, comprehensive. Besides such topics as 
the distribution, growth, forms, and physical, chemical, and mechani- 
cal properties of sea ice, he dealt in detail with the ocean/atmosphere 
system and its relationship to the pack ice. In addition to qualita- 
tive description of the physical events, he often gave them a quanti- 
tative formulation through the use of empirical or carefully derived 
mathematical expressions. Although subsequent research has shown some 
of Zubov' s work to contain oversimplifications, many of his contribu- 
tions remain of value and his text is still regarded as a standard 
reference. 

Emphasis during early arctic reasearch was upon observing and 
describing events occurring on the ice and in the ocean. Not until 
the early 1950 's did the concept of relating the energy balance to 
the mass budget begin to gain favor. Earlier there had been only 
sporadic attempts to measure the radiative components, but, with the 
establishment of NP-II in 1950 and T-3 in 1952, began aserious effort 
to acquire representative data on both short-wave and long-wave radia- 
tion. These measurements , combined with calculations of the turbulent 



heat fluxes, produced for the first time good quantitative estimates 
of the heat-transfer processes in the Central Arctic. The results of 
a decade of observations have been summarized by Gavrilova (1963). 
Not only has she presented the geographical distribution of the radia- 
tion components in the Arctic and peripheral areas, but she has also 
described Soviet techniques of observation and methods of analysis. 
Because the Soviet data are supported by the most field observations 
and thus seem to be the most representative, Gavrilova's summary must 
be considered as a prime source of data concerning the arctic radia- 
tion regime. 

In 1962 the Arctic Meteorology Research Group at McGill Univer- 
sity, led by S. Orvig, undertook a project to correlate and evaluate 
all available data relating to the arctic heat budget. Included in 
these studies were such topics as atmospheric advection of heat and 
moisture, heat flux through the ice, sensible and latent heat fluxes, 
radiative components, and surface albedos. Despite certain limita- 
tions (see Section 4.1.4), the resulting technical papers, published 
during the next four years, represent a significant contribution to 
knowledge of the arctic regions. 

Pursuing a similar theme, Fletcher (1965) related the arctic heat 
budget to hemispheric climate. In this study he discusses the inter- 
action of the arctic ice mass with the ocean and atmosphere, mentions 
possible methods by which man may influence the heat balance, comments 
on the present direction of arctic research, and specifies the re- 
search topics deserving priority. One aim of this study was to obtain 
an average heat budget for the central Arctic Basin today for compari- 
son with a possible ice-free regime. After a detailed analysis, he 
proposed a heat budget that relied heavily on Soviet sources, but 
which did not correspond exactly with any other single author. This 
budget is internally consistent and, according to Fletcher, has a pro- 
bable error less than the year-to-year variations. Much of the work 
described in the present report makes use of these values. By inter- 
polating cloud and energy exchange data from the normally ice-free 
portions of arctic seas, Fletcher devises a probable energy budget 
for an ice-free central Arctic Basin. The net energy balance for this 
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situation is very slightly positive, suggesting that the ice would not 
immediately re-establish itself once it is removed. However, result- 
ing changes in atmospheric circulation could tip this delicate balance 
in either direction. Thus the stability of an ice-free Arctic remains 
an open question. 

1.5 UNSOLVED PROBLEMS 

It has long been observed that unusual weather in Europe seems 
to be correlated with abnormal extent of the pack ice. This feature 
is to be expected since storms generally tend to follow the margins 
of the ice, drawing energy from the temperature difference between 
ice and water. In an attempt to define the factors influencing cli- 
matic change, for example, Lamb (1959, 1961) made a detailed study of 
world-wide weather records and found that during the period 1800 to 
1940 there was a definite tendency toward a stronger general circula- 
tion. The edge of the pack, ice at the same time showed such a large 
northward migration that by 1938 it had receded further north than 
ever before in modern times. Whether the abnormally small amount of 
sea ice caused the circulation to be more vigorous, or whether the 
greater heat transport by the atmosphere was responsible for less ice, 
is a question that has not been adequately answered. Lamb attributes 
the recession of the ice to increased atmospheric heat advection as- 
sociated with a more vigorous global circulation. More recently, 
Fletcher (1969) has suggested that the fluctuations in global atmos- 
pheric circulation and in the state of the ice and general circulation 
in the Northern Hemisphere may be related to the great changes in ice 
cover which occur in the Antarctic, but data remain so sparse that 
evaluation is very difficult. Until theoretical models describing 
the entire ocean/ice/atmosphere system can be solved, all that may be 
definitely stated is that some causal relationship must exist between 
the ice extent and hemispheric weather patterns. 

Since relatively small changes in the energy balance over the 
pack ice can produce large changes in the areal extent of the ice 
cover, numerous investigators have looked to the Arctic for a source 
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of the sweeping climatic changes associated with the ice ages of the 
Quaternary. The present pattern of atmospheric and oceanic circula- 
tion and the distribution of ice may be only one of several quasi- 
stable situations. This possibility has evoked speculation on the 
various consequences of an ice-free Arctic Ocean. Many investigators 
conclude that an open Arctic, if it indeed represents an equilibrium 
situation, could lead to the accumulation of continental ice sheets. 
The central question to be answered is then: "What natural processes 
could produce an ice-free Arctic, and, once an open Arctic is achieved, 
could it remain so?" Brooks (1949) postulated a critical size of the 
ice pack above which its own cooling effect would ensure its continued 
growth, but, below which, a warming climate would ensure its continued 
wastage. Budyko (1966) presents a similar positive feedback mechanism 
to explain the formation of an open Arctic Ocean. In their theory of 
the ice ages, Ewing and Donn (1956, 1958, 1966) advanced data from a 
multitude of sources arguing for the existence of an ice-free Arctic, 
concurrent with the continental ice sheets. They suggest that varia- 
tions in the oceanic heat flux, caused by a rising and falling sea 
level, could explain the appearance and disappearance of the pack ice. 
Because it is not known how large the oceanic heat flux must become 
before the ice would vanish, it is not clear whether the rise in sea 
level, associated with the melting of continental glaciers, could in- 
crease the supply of energy enough to melt the ice. In fact, recent 
evidence (Hunkins , 1965; Ku, 1965) strongly indicates that sedimenta- 
tion rates on the floor of the Arctic Basin have remained relatively 
constant over at least the last 70,000 years. If so, it would appear 
that the Wisconsin Glaciation occurred without a coincident open Arc- 
tic. 

As man's technological capabilities increase, it becomes more 
and more urgent that the factors influencing climate on a large scale 
be understood. For example, it now may be possible for man to remove 
the arctic ice pack; the result can only be surmised. Another possi- 
bility is that in man's increasingly effective tampering with nature, 
he may accidentally intervene at some critical stage in the climatic 
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process, producing unexpected effects on a world-wide scale. Such 
considerations cannot be neglected in the formulation of climatic ex- 
periments . 

Because the physical principles governing the operation of each 
part of the system are relatively well understood, theoretical models 
of each part can be designed and methods of solution devised. The 
complexity of the resulting equations usually necessitates a numeri- 
cal solution, employing the use of high-speed computers. Although 
storage limitations in present-day computers make it possible to solve 
only simplified versions of each component of the system, forthcoming 
computers offer hope that more realistic models may be solved. The 
ultimate goal then must be the incorporation of all the individual 
models into a single program that will provide a simultaneous descrip- 
tion of each component. 
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II. FORMULATION OF A THEORETICAL MODEL 

2.1 APPROACH TO THE PROBLEM 

Among the many unsolved problems concerning the Arctic, that of 
quantitatively predicting the response of sea ice to changes in its 
environment stands out. Not only is it of major importance, but it 
also is a relatively straightforward problem amenable to current com- 
puter techniques. However, since theoretical treatment of the atmos- 
phere and the ocean is beyond the scope of this report, their effects 
must be specified as external parameters. The outcome will be that 
the environment may influence the ice, but not vice versa. The ob- 
jectives of this study will then be to forecast ice temperatures and 
thickness and to define the role of each component of the energy bud- 
get in relation to its effects upon the ice. 

2.1.1 Previous work 

A large number of studies have been made to develop empirical 
relationships predicting ice growth from observed temperatures at the 
surface. Two of the best known empirical formulas are those of Barnes 
(1928) and Zubov (1938). Somewhat more complex equations have been 
developed by Bilello (1961, 1964) from long-term observations of ice 
growth in the Arctic. While empirical and statistical approaches have 
practical value, they offer little physical insight into the problem. 

Since the classical work of Stefan (1891) , theoretical studies of 
the rate of ice growth have generally been limited to those in which 
simple analytical, approximate solutions could be attained. One of 
the most successful of these is due to Kolesnikov (1946) . His equation 
predicts the thickness of sea ice by setting up heat-balance equations 
and solving them in conjunction with the heat conduction equation. 
Solution requires specification of such parameters as wind speed, hu- 
midity, air temperature, cloudiness, snow thickness, surface roughness, 
and ice salinity. This equation predicts growth rates in almost exact 
agreement with observed ones ( Calloway , 1954). In practice, it would 
be difficult to acquire values for all these parameters over a large 
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region. Doronin (1959, 1966) calculated growth rates for an ice slab 
by using the heat-conduction equation, allowing for the variation of 
specific heat with ice salinity. He reports good agreement with ob- 
servations on thin ice, but gives no specific examples. 

In a more sophisticated model, Kolesnikov (1958) proposed a set 
of steady-state Fourier equations to describe the air/snow/ice/water 
system; approximate analytical solutions were then given. However, 
several serious errors have been pointed out in the model (U. S. Naval 
Civil Engineering Lab . , 1965), and hence its applicability is doubtful. 

With the development of high-speed computers, it is no longer 
mandatory that models have analytical solutions. A system of integro- 
differential equations can now be rapidly solved by finite-difference 
techniques. An advanced model is reported by Budyko (1966). Using 
a method of successive approximations, he reports that he has calcu- 
lated ice thicknesses and the horizontal extent of the ice pack in 
good agreement with observational evidence. This model requires speci- 
fying the heat flux from the ocean at the bottom of the ice; abla- 
tion at the upper surface is determined by assuming values for humid- 
ity, wind speed, ice temperatures, air temperatures, and the radiative 
fluxes. The principal criticism of Budyko 's model is that its upper 
boundary condition is more or less physically independent of heat con- 
duction within the ice. By specifying surface temperatures, Budyko 
ignores any effects of heat conduction on surface ablation and surface 
temperature. No provisions are made to account for snow cover or the 
effects of ice salinity on heat transport between the boundaries of 
the slab. A model described by Untersteiner (1966) eliminates all 
these difficulties. 

2.1.2 Idealized picture of the ice 

In an idealized fashion, arctic sea ice may be viewed as an in- 
finite, horizontally homogeneous slab, floating on its own liquid 
phase, having various energy fluxes impinging at both surfaces. The 
ice transfers heat between ocean and atmosphere by conduction; however, 
the transfer is affected by the brine pockets in the ice and by short- 
wave radiation penetrating the upper surface during certain seasons. 
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The annual snow cover influences the radiative and conductive 
transfer in important ways. It is taken into account by addition of 
an appropriate second layer to the model. An internal boundary con- 
dition ties together the diffusion equations describing conduction of 
heat in the snow and in the ice. 

The boundaries of the slab are pictured as mathematically ideal- 
ized planes. Energy is assumed to be absorbed at these planes, al- 
though in reality the energy absorption occurs in a finite layer. 
Mass changes may occur at either the upper or the lower surface. At 
the upper boundary, ice or snow may melt, but mass may be increased 
only through addition of snow. Either ablation or accretion can take 
place at the bottom of the slab. Whether ablation or accretion occurs 
is determined by the energy balance at the boundary under considera- 
tion. If, for example, there is more energy being absorbed at the 
boundary than is being lost from it, then either ablation or warming 
must occur. 

The external energy fluxes, which are basic in determining mass 
changes of the slab, are taken to repeat themselves in a constant an- 
nual cycle. Thus an equilibrium situation must eventually be reached, 
that is, the net yearly bottom accretion will exactly balance the 
amount of ice lost through melting at the upper surface. This means 
that the net mass balance for the year is zero and that the mass 
changes and temperature patterns will begin to exactly repeat them- 
selves in an annual cycle. Either condition may be used to define 
equilibrium. 

Conditions are often uniform over vast areas in the Arctic. 
Available radiation data (Gavrilova , 1953) show only a slight latitu- 
dinal variation over the ice-covered portions of the ocean. Thus the 
assumption of horizontal homogeneity appears to be justified. 
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2.2 HEAT CONDUCTION 

2.2.1 Heat conduction in the ice 

The basic equation governing heat flow in solids is the thermal 
diffusion equation: 

PC ~ = V • (kVT), (1) 

where p is the density, c the specific heat, T the temperature, and 
k the thermal conductivity, and t refers to time. This well known 
equation has received intensive study, and numerous analytical solu- 
tions have been developed for certain situations ( Carslaw and Jaeger , 
1959). However, because of the deep penetration and absorption of 
short-wave radiation, Eq. (1) does not completely describe the situa- 
tion as it exists for sea ice. Untersteiner (1964) used the following 
modifications of Eq. (1): 

(pc) i f = h < k i f > + "±* 0 e *P <-i z > - (P'>i * £ , <*> 

where is the extinction coefficient, I q is the amount of short-wave 
radiation passing through the surface, and w is the upward flux of ice 
due to thickness changes. The subscript "i" refers to variables in 
the ice. The third term in Eq . (2) represents the heat generated by 
absorption of solar radiation per unit of time per unit volume at a 
depth z, and the last term describes the heat flux associated with the 
hydrostatic readjustment of the ice in response to mass changes at the 
boundaries. In keeping with the assumption of horizontal uniformity, 
the problem was reduced to one spatial dimension. 

Both (pc)_^ and k^ are variables in sea ice, because small pockets 
of brine, trapped during freezing, remain in the ice. The brine is 
assumed to be at its freezing point and in phase equilibrium with the 
surrounding ice. The equilibrium is maintained by volume changes in 
the brine pocket. A rise in temperature causes the ice surrounding 
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the pocket to melt, diluting the brine and raising its freezing point 
to the new temperature. The brine pocket, thus, is a thermal reser- 
voir, retarding the heating or the cooling of the ice. It is obvious, 
hence, that density, specific heat, and thermal conductivity are all 
functions not only of salinity, but also of temperature. 

The variations of (pc)^ and k^ were first investigated by 
Malmgren (1927), and a thorough theoretical analysis of the problem 
has been presented by Schwerdtf eger (1963). Untersteiner (1961) in- 
troduced the following formulas to describe these variations: 

fncV = focV .+ T s(z) _ m 



where S(z) is the ice salinity at a depth z; y « 4100 cal °K/gm and 
g = 0.28 cal cm /gm sec are constants. The subscript "f" refers to 
pure ice and (pc^ f - 0.45 cal cm 3 °K and ^ f = 0.00486 cal/cm sec °K. 
Ono (1966) derived expressions for (pc)^ and k^ on theoretical grounds. 
The resulting equations, while slightly more complex than Untersteiner ' s 
proved to be similar in form and to give values which are in remarkably 
good agreement. For simplicity, Eqs . (3) and (4) were chosen for use 
in this model. 

Substituting Eqs. (3) and (4) into (2), Untersteiner (1964) de- 
rived a complete description of energy transfer within the ice: 



Y S(z) 



3T 



. bs(z) ri. 

T - 273 .21 



dS 3T 



(PC), 



(5) 
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Solution of this equation, with specified boundary temperatures, pro- 
vided a temperature field consistent with observations. 

For the present experiment, Eq. (5) is cumbersome and its solu- 
tion time consuming, but an examination of the equations reveals that 
several terms may be dropped without substantial loss of accuracy. 
The last term in Eq. (5), the vertical advection term, is small and 
was eliminated. Its removal has the effect of allowing the upper and 
lower boundaries to move independently of one another; i.e., there is 
no hydrostatic readjustment. The terms - ,^273 |~ and 

6 S(z) 9T 2 ^ . . 

9 express the variatxon of conductivity within the layer 

(t - 273r dz 

under consideration. Because these terms are awkward to handle, we 
decided to assume that the conductivity was constant within an incre- 
mental layer of thickness Az. This assumption still allows the equa- 
tion to take into account the variation of conductivity from layer to 
layer. An order-of -magnitude analysis performed on the equation indi- 
cated that both terms are three orders of magnitude smaller than the 
major terms; hence their removal seemed justified. An experiment de- 
scribed in the next chapter further serves to demonstrate that their 
effect is truly negligible. 

With the above simplifications, Eq. (5) becomes: 

L c <v + I s ( z > 1 il _ L . g S(z) 1 S^T 



Eq. (6) is now taken to be the basic equation describing hea 
in the ice. 



2.2.2 Heat conduction in the snow 

An analogous equation may be employed to describe passage of heat 
through a snow cover. Assuming horizontal homogeneity and a constant 
thermal conductivity, Eq. (1) becomes for the snow: 
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(7) 



where the subscript "s" refers to values in the snow. 

Snow, unlike ice, undergoes great changes in density and crystal 
fabric during its metamorphism. Also the extinction coefficient for 
short-wave radiation varies widely. Observations by Liljequist (1956) 
and Ambach (1962) indicate that < s may range from 0.07 cm" 1 to 0.5 cm -1 , 
depending on temperature, density ^ and crystal structure. A labora- 
tory experiment by Mellor (1965) produced values between 1.3 cm -1 and 
1.7 cm for cold snow. The consequence of such large values is that 
essentially all incoming radiation is absorbed in the first few centi- 
meters of the snow. Thus, in the model, the heating effect from this 
source is indistinguishable from absorption at the surface. In the 

absence of any detailed information on the behavior of < , heat con- 

s 

duction in the snow may be described simply by: 



2.3 BOUNDARY CONDITIONS 

The crux of the present problem is the description of mass changes 
at the boundaries of the slab. They are stated as heat-balance equa- 
tions , relating the various energy fluxes toward and away from the sur- 
face to the rate of freezing or melting. Solution of the boundary 
equations, along with the two heat-conduction equations, will allow 
computation of surface temperatures, temperature profiles within the 
ice and snow layers, and mass changes at both boundaries of the slab. 
The present model thus will represent a generalization of Untersteiner ' s 
(1964) model. 

2.3.1 Upper surface 

The major energy fluxes at the upper boundary include incoming 
long-wave radiation from the atmosphere and clouds (F ) , incoming 




(8) 
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short-wave radiation (F^) > reflected short-wave radiation from the 
surface (aF^) , and outgoing long-wave radiation (e^ • * n Edi- 

tion to these, there are several smaller but important fluxes: the 
fluxes of sensible heat (F g ) and latent heat (F^) , heat conduction 
flux in the ice or snow ( F c ) » and a flux of radiative energy through 
the surface into the ice (I Q ) • These fluxes are schematically illus- 
trated in Fig. 2. In describing the energy fluxes, the convention 
will be adopted that a flux toward the surface is positive and one 
avay from the surface is negative. 

To determine a balance equation for the upper boundary, it is 
necessary to consider two possible situations. If the surface tem- 
perature (T q ) is below the freezing point, then T q must be adjusted 

in order that all the fluxes balance. If, on the other hand, T is 

o 

at the freezing point, a certain amount of ice will have to melt to 
accommodate any surplus of energy flux toward the surface: 



where h is the thickness of the snow layer, H is the thickness of the 
ice layer, is the long-wave emissivity, a is the Stef an-Boltzmann 
constant, q is the latent heat of fusion, and a is the surface albedo 



of the ice, it represents an energy loss by the surface and must be 
treated accordingly (see Section 2.4.2). The subscript "o" refers to 
the upper surface. The turbulent fluxes have been given as positive 
terms, but they may in fact be in either direction. The quantities 
in Eq. (9) which must be specified as time-dependent external para- 
meters are a, F^, F^, I q , F g , and F^; all temperatures and ice-thick- 
ness changes result from integration of the model. Snow deposition 
must also be treated as an external parameter (see Section 2. 4. A). 




(9) 



I q is contained in (1-cOF^; hi 



, since I passes into the interior 
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Although Eq. (9) is quite general, it is only one of several pos- 
sible formulations. For example, instead of specifying the downward 
long-wave and short-wave radiation at the surface, one might develop 
the equation in terms of solar radiation incident at the top of the 
atmosphere and assumed conditions of cloudiness, atmospheric absorp- 
tion, humidity, and temperature. Added complexities of this sort, how- 
ever, were not felt to be justified at the present stage of develop- 
ment. 

2.3.2 Ice/water interface 

At the bottom of the ice there are only two fluxes to be consid- 
ered: the turbulent heat flux from the ocean and the conductive heat 
flux in the ice close to the boundary. The magnitude of the two fluxes 
is similar, with one or the other dominating during different seasons. 
The boundary equation has the form: 

3T 3T , 

*i IT h+H " K » >T h+H - 'It <h+ » h+H • <«» 

where K is the coefficient of eddy diffusivity in the water, and the 
subscript "w" refers to values in the water; h and H are the thick- 
nesses of snow and ice. Without a theoretical model of the ocean, it 
is not possible to calculate the heat flux from the water. Thus 

(pc) w (K w 

growth is determined by the model. As a convenience, the following 
notation will be adopted throughout the remainder of the text: 
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2.3.3 Ice/snow interface 

When snow covers the ice, it becomes necessary to specify another 
restriction at the ice/snow interface (z = h) : 



i 8z h s 3z h 

Equation (12) simply states that conduction through the ice/snow inter- 
face is continuous and assumes that any effects from penetration of 
short-wave radiation are negligible down to the depth h. Also, it is 
required that at z = h, T. = T g . 

In summary, the problem is to solve the four equations (6), (8), 
(9), and (10) for the four unknowns: snow temperature, ice tempera- 
ture, thickness change at the top, and thickness change at the bottom. 
This model, which is equivalent to the one outlined by Untersteiner 
(1966), is presented graphically in Fig. 2. 

2.4 EXTERNAL PARAMETERS 

The model formulated in the preceeding sections views temperature 
and thickness changes as a result of arbitrarily imposed energy fluxes 
at the upper and lower boundaries of the slab. Solution of the equa- 
tions depends upon specification of these fluxes, along with a number 
of other parameters. In the initial integration of the model, values 
which closely resemble actual conditions were chosen; results could 
then be compared with observations as a test of the validity of the 
model. 

2.4.1 Energy fluxes 

At the upper boundary, four energy fluxes must be specified: 
F r' F L' F s' and F l' After a detailed analysis, Fletcher (1965) chose 
the values given by Marsh un ova (1961) for and F^ as the most repre- 
sentative of the region. Her values were based on the results from 
Soviet drifting stations NP-3 and NP-6. Fletcher also considers the 
turbulent fluxes calculated by Doronin (1963) physically plausible and 
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consistent with what is known about the other components. Following 
Fletcher's suggestion, these values were adopted for use in the ini- 
tial test of the model and are reproduced in Table 1. The only major 
researchers to deviate substantially from the above findings are 
Vowinckel and Orvig (1966, 1967). Their studies indicate values for 
the absorbed short-wave radiation which seem extremely high. The con- 
sequences of the heat budget of Vowinckel and Orvig are examined in 
detail in Section 4.1.4. 

Compared with the fluxes at the upper surface, the magnitude of 
the oceanic heat flux (F ) is relatively little known. Direct measure- 
ments have proved difficult, at best, and estimates have been based 
primarily on indirect evidence. Malmgren (1927) used temperature pro- 
files in the ice to produce his estimate of 6.8 kcal/cm 2 year. How- 
ever, equating the annual heat flux in the ice with that in the water 
is incorrect (Untersteiner , 1964) , since this method ignores the heat 
released from accretion at the underside of the ice and heat contri- 
buted by the penetration of visible radiation (similarly Mosby , 1963) . 

One method of approximating F w is to estimate the heat content 

and the residence time of the waters entering the Arctic Basin. Rough 

values of the heat loss can then be calculated for various areas. 

Crary (1960) used the data from Timofeev (1958) to construct lines of 

equal heat content for the North Polar Basin. His chart indicates a 
2 

heat loss of approximately 2 kcal/cm year from the ocean in the Cen- 
tral Arctic. A more recent calculation (Panov, 1964) demonstrates 
that the heat flux undergoes large year-to-year variations. In the 
Chukchi Sea, for example, the annual upward heat loss from Atlantic 

waters fluctuated between 1038 and 1800 -10 12 kcal during the 1930 to 1964 
12 

period; the average being 1426 -10 kcal/year. Assuming an area of 

2 2 
86,400 km , the annual heat flux then varied from 1.1 to 2.0 kcal/cm , 

with an average of 1.65 kcal/cm^. It would be reasonable to expect a 

slightly greater value in the Central Arctic. 

Badgley (1961) used a more direct method to determine F . By 

choosing a period when the lower part of the ice was isothermal and 

then measuring bottom ablation, he calculated a yearly value of 1.5 

kcal/cm 2 . Independent results from Untersteiner 's 1964 model suggest 
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almost exactly the same value. Although all these methods involve 
numerous assumptions, the true values of F w must lie close to 1.5 
kcal/cm^ year. This value was adopted for the test run of the model. 

A final problem connected with heat flux from the ocean is that 
there is probably a slight time dependence. During the summer, fresh 
water drains into the ocean and forms a stable layer near the under- 
side of the ice (Unters teiner and Badgley , 1958; Hanson , 1965). But 
without more precise knowledge than is currently available, it seems 
pointless to speculate on time variations of 7^. Therefore it was 
assumed to be constant during the entire year. 

2.4.2 Penetration by short-wave radiation 

During the snow- free summer month's , solar radiation penetrating 
the ice surface (I q ) is the prime source of internal heating. In this 
section we will be concerned with the amount of penetrating radiation 
and its effect. 

Temperature profiles taken at the height of the melt season in- 
dicate an internal heating of the slab of approximately 700 cal/cm 2 
in only 14 days ( Unters teiner , 1961) . A further examination of these 
profiles shows that the temperature gradients (and therefore the con- 
ductive heat fluxes) remained nearly constant during this period, im- 
plying that the observed heating was not the result of conduction. 
Thus the source of the energy used to heat the deeper layers must be 
ascribed to penetration of incoming short-wave radiation. Unter- 
steiner's data indicate that about 32 percent of the net short-wave 
radiation makes no contribution to the surface melting. On the basis 
of these observations and his radiation measurements , Unters teiner 
(1964) arrived at an annual value for I of 1.2 kcal/cm . 

Because absorption of solar energy by the ice is a function of 
wavelength, angle of incidence, and physical character of the surface, 
a theoretical description of I q would be complex. Ice is nearly 
opaque to some wavelengths , and the incoming energy at these wave- 
lengths is absorbed in the top few centimeters. At other wavelengths, 
the radiation penetrates farther. Thus, for modeling purposes, we 
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may consider the ice slab as having two layers: a thin surface layer 
that absorbs most of the short-wave energy, overlying another layer 
in which the remaining energy is absorbed. Numerical limitations de- 
crease the resolution of the model to the point that absorption of 
radiation in this surface layer is indistinguishable from direct sur- 
face melting. As a first approximation, we may assume that Beer's 
Law is valid within the lower layer and determine a bulk (integrated 
over all wavelengths and solid angles) extinction coefficient (k^) . 
Direct photometric measurements ( Chernigovskii , 1963) indicate values 
for K of 0.014 cm 1 to a depth of 150 cm, and 0.012 cm" 1 below this. 
Untersteiner (1961) utilized temperature profiles to make indirect 
calculations of energy transmission in the ice. His results give an 
average value of 0.015 cm -1 between 50 and 150 cm. Since both of 
these methods have shown that extinction in the lower layer closely 
follows Beer's Law, Eqs . (6) and (9) offer a reasonable description 
of the effects of penetration and absorption of short-wave radiation. 

Although could easily be specified as a function of depth in 
the model, little is known of such variations and a constant value of 
0.015 cm 1 was therefore chosen. To specify I q in the model, values 
for incoming short-wave radiation were combined with the albedo to 
calculate net short-wave radiation; a predetermined percentage of 
this quantity was then assigned to I q (subtracted from the net short- 
wave radiation) and the remaining portion was used to determine the 
surface energy balance. Lacking more detailed evidence, this percent- 
age was taken to be constant during the snow-free period. Imposing 
the restriction that I q be 1.2 kcal/cm 2 year on Fletcher's suggested 
heat budget means that only 17 percent of the net short-wave radia- 
tion may penetrate the ice. This is in apparent conflict with Unter- 
steiner *s 32 percent and it is possible that his 1964 value under- 
estimates I q . The model's predictions of surface ablation and ice 
thickness should help to resolve these discrepancies. 
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Table 1. AVERAGE MONTHLY VALUES FOR THE INPUT DATA AT THE UPPER BOUNDARY 
i ACCORDING TO FLETCHER (1965) . 
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0.83 0.81 0.82 0.78 
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2.4.3 Ice salinity 

Brine inclusions in sea ice are thermodynamically important be- 
cause they affect density, conductivity, and latent heat. Equations 
(3) and (4) describe the influence of salinity on (pc^ and and 
hence upon heat conduction. At temperatures near the freezing point, 
(pc)^ may change by almost two orders of magnitude, while shows a 
10-percent change. Also the latent heat of fusion (q) is influenced 
by the initial brine volume. At the bottom of the ice, for example, 
the brine volume is roughly 10 percent ( Schwarzacher , 1959) , and the 
heat of fusion 64 cal/cm (rather than 72 cal/cm 3 ) . This difference 
is important in the calculation of bottom growth. 

Although it is relatively easy to treat the effects of salinity, 
when it is known, acquisition of representative salinity profiles has 
proved difficult because of the large scatter in the field data. 
To further compound the problem, the shapes of the profiles probably 
depend upon the season, the growth rate, and the thermal history of 
the slab. Many of these difficulties could be resolved by a theoret- 
ical description of the desalination process, however at the present 
time, there is still doubt as to what mechanism controls brine re- 
moval. A recent study ( Untersteiner , 1968) has shown that solute 
diffusion and gravity drainage are too small to account for the ob- 
served features. Both "flushing" and "brine expulsion" could produce 
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br ine withdrawal at a reasonable rate, but it is now known which is 
dominant . 

Ideally, salinity should be described as a function of growth 
rate and thermal history of the slab; however, a proven theory of de- 
salination is lacking, and it is difficult to guess how the profile 
would change under varying conditions of temperature and speed of 
growth. Fortunately an examination of Eqs. (3) and (4) indicates 
that unless the temperature is higher than -2 or -3°C, salinity has 
only a small effect upon (pc) ± and V Thus it is important only that 
salinity be known accurately during the summer and in the vicinity of 
the ice/water interface. Salinity in the lower portion of the ice is 
relatively constant at around 3.2 parts per thousand. The best doc- 
umented equilibrium profile is that of Schw^acher (1959) shown as 
Fig . 3. In the absence of more extensive data, this profile was adop- 
ted to represent mean annual conditions . 
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Fig. 3 -- Salinity profile for equilibrium sea ice. 



The salinity of Fig. 3 is approximated by the expression: 



t) = A + B sin 



, n/(f + m) 



where A, B, C, D, m, and n are constants and H is the thickness of the 
ice. Equation (13) thus gives the salinity as a function of depth and 
is related to time through H. Details of the evaluation of the con- 
stants are given in the following chapter. The equation, while com- 
plex, provided an extremely good fit to the data. Implicit in the 
equation is the assumption that shape and end-points of the profile 
are retained during all fluctuations of ice thickness. 

2.4.4 Snow cover 

Nearly two decades of observations from various drifting stations 
have produced a fairly good knowledge of the snow distribution in time 
and space. After the freeze-up in late August, sriow accumulates rap- 
idly for the first month or two, followed by a more gradual increase 
during the winter months. In May, snow accumulation again slightly 
increases , reaching a maximum depth of about 40 cm in June ( Unter- 
steiner , 1961; Hanson , 1965). During the melt season, sporadic storms 
deposit small amounts of snow that directly influence the mass budget 
only slightly, but which, by changing the albedo, may affect the heat 
and mass budget considerably (see Section 2.4.5). Although surface 
features, such as pressure ridges, cause some small-scale variations, 
snow depth is generally uniform throughout the Central Arctic. 

To represent this pattern of snowfall, it was assumed that there 
is a linear accumulation of 30 cm between 20 August and 30 October, 
a linear increase of 5 cm from 1 November to 30 April, and an addi- 
tional 5 cm during the month of May. In cases where the summer melt 
season extends past 20 August, snow accumulation is delayed until 
freeze-up occurs. If the surface temperature reaches the melting 
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point before 1 June, then snow accumulation is cut off before it 
reaches 40 cm. No attempt was made to model the occasional summer 
snowfalls . 

The fine-grained structure of arctic snow, combined with com- 
paction by the wind, causes it to have a fairly high density. A value 
3 

of 0.33 g/cm ( Untersteiner , 1961) was chosen for the model. It is 
assumed that the snow density (p g ) is independent of depth and remains 
constant as long as the snow temperature remains below freezing. At 
the onset of melting, the density increases rapidly to a limiting val- 
ue of 0.45 g/cm (E. LaChappelle, personal communication) as melt 
water penetrates the snow and is refrozen. Thermal conductivity (k ) 
is quite variable, even in snows of equal density, and depends upon 
grain size, snow structure, and the temperature gradient. Mellor 
(1964) has collected and graphed experimental results from 75 years 
of observations. Abels' Formula ( Abels , 1892) represents a reason- 
able compromise of all the diverse measurements and was used in the 
experiment : 

k g = 0.0068 p^(cal/cm sec °K) . (14) 

A major problem encountered in treating the snow cover is how to 
take into account water percolation and densif ication of the melting 
snow. It would be possible, for example, to add a forcing function 
term to Eq. (8) that would describe the heat released by refreezing. 
Equation (14) could then be combined with an assumed p g to solve the 
modified heat-conduction equation. Before this method was tried, a 
calculation was made to determine the heat content of the snow layer 
at the onset of melting. The results indicated that the melting and 
subsequent refreezing of only 2 cm of snow would add enough heat to 
make the entire layer isothermal. No field data are available on the 
length of time necessary to melt the first 2 cm of snow; however 
an initial test of the model predicted that only about 5 days are 
required. In order to avoid unnecessary complications, the model 
was then set up to ignore heat added by melt water. The integration 
proceeds, considering only conduction in the snow, until 2 cm 
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of snow are melted. At this time, the snow layer is made isothermal 
and the density is increased from 0.33 to 0.45, with a consequent dis- 
continuity in snow depth. This procedure preserves the essential fea- 
tures of the densification process. The only significant deviation 
from reality is that the temperature at the ice/snow boundary will be 
too cold for a few days and then will suddenly jump to the freezing 
point of the ice (-0.1°C). Consequences of this treatment are dis- 
cussed more thoroughly in Section 2.5. 

2.4.5 Surface albedo 

The surface albedo (a) and its variations are probably the most 
important regional factor affecting the heat and mass budgets of the 
Arctic Basin. Albedo values determine the net short-wave radiation, 
and small changes in a can account for dramatic changes in the amount 
of surface ablation. As an example, if the albedo were decreased by 
only 10 percent during the month of June, yearly ice ablation would 
increase by about 25 cm or 60 percent. 

During the spring and fall, when surface conditions are uniform, 
a is high and relatively constant. Incident short-wave radiation is 
relatively small during this time. The critical period occurs during 
June and July, when there is a sharp increase in available energy and 
melting is rapid. As melting progresses, wide albedo fluctuations 
occur in response to the changing micro-relief of the surface. Yakov- 
lev (1958) reports that a varied between 0.44 and 0.89 on NP 2 during 
the summer ablation season. 

Snow patches, bare ice, and melt ponds, all with different radia- 
tive characteristics, may coexist over a relatively small area. In 
addition, snowfall can cause a 30-percent increase in a in a matter 
of hours. It is obviously impossible to duplicate such random fluc- 
tuations and horizontal variations in a one-dimensional model. 

Most authors agree on spring and fall values for a, but there is 
disagreement regarding albedo averages during the summer. The major- 
ity view is that the effective albedo of ice that is melting lies be- 
tween 0.60 and 0.70. Untersteiner (1961) assumed an average of 0.66, 
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while Soviet measurements indicate a value of 0.64 (Marsh unova , 1961; 
Chernigovskii , 1963) . Aircraft measurements ( K. Hanson , 1963) also 
support these values . Unpublished data from Arlis-II show hourly al- 
bedo values fluctuating from 0.55 to 0.70 on melting ice. 

Vowinckel and Orvig (1966, 1967) hold a value of 0.45 to be rep- 
resentative. This view is supported by the fact that albedo measure- 
ments on melting, dirty ice have often been as low as 0.40 and inte- 
grated albedo values from a 15 meter tower are even lower ( Langleben , 
i968). These facts, combined with the prevalence of melt ponds, cause 
Vowinckel 's estimate to appear reasonable: however acceptance of a 
value this low is incompatible with other considerations. From energy- 
balance calculations, the increase in net short-wave radiation result- 
ing from an albedo of 0.45 would cause about 120 cm of ablation during 
the summer, an amount roughly three times the observed quantity. On 
this basis, Fletcher (1965) accepted Marshunova's value of 0.64. 

Neglecting the influence of melt ponds, the average ice albedo 
is certainly closer to 0.64 than to 0.45, although an areal average 
is probably substantially lower. The question which must be consid- 
ered is how the depth of the melt ponds affects the average ice thick- 
ness over an annual cycle. A melt pond which exists at the end of the 
ablation season will freeze during the fall and, except for acting as 
a heat source to retard cooling, has little effect on annual thick- 
ness. During the summer, however, at least half of the melt ponds 
drain, forming channels and hollows in the surface topography. Snow 
accumulated in these depressions remains longer, retards ablation in 
these areas, and tends to level the ice relief. Increased ablation 
beneath melt ponds therefore has only a limited effect on the general 
ice thickness. In consequence, averages over large areas were set 
aside, and the albedo of melting ice was taken as 0.64. Fall and spring 
snow albedos were also taken from Marshunova (1961) and are reproduced 
in Table 1. 

A further problem is how to model the albedo during the transi- 
tion from a snow cover to bare ice. Lacking any detailed information, 
a simple linear decay of albedo with snow depth was assumed: 
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0 . 0i + (8o - a±) (15) 

where h is the snow depth when the snow begins to melt, a Q is the 
snow albedo at the onset of melting, and a ± is the albedo of the bare 
ice (0.64). Thus, if the snow albedo was 0.84 at the onset of melting, 
it would have a value of 0.74 when half the snow had melted, 0.69 after 
three-quarters had melted, etc. Although it is not known how well 
Eq. (15) approximates reality, the errors introduced should be no 
larger than those from other sources. 

2.4.6 Other assumptions 

For convenience, the year was divided into twelve months of 30 
days each. Since both ice and snow radiate in the long-wave region 
much like a black body ( Kellogg , 1964) , the long-wave emissivity (e^) 
in Eq. (9) was set equal to unity. It may be seen from Eqs. (3) and 
(4) that both (pc) ± and k^, are discontinuous at T = 273°K. In order 
to avoid this singularity, the melting temperature was assumed to be" 
272. 9°K, which is realistic, as a residue of salt remains in the upper 
layers of the ice. The temperature at the underside of the ice was 
assumed to remain constant at the freezing point of sea water (271. 2°K). 

2.5 LIMITATIONS 

There are a number of limitations that are inherent to, or implicit 
in, the proposed model. These shortcomings may be roughly divided into 
two groups: conceptual and physical. Conceptual limitations arise from 
a lack of theoretical understanding or the inability to model certain 
processes affecting the ice; physical limitations are due to uncertain- 
ties in the environmental data. The following discussion will outline 
these difficulties and attempt to evaluate the importance of each to the 
purpose of the present calculations. 
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2.5.1 Conceptual limitations 

A basic shortcoming of the model is its inability to account for 
effects of mechanical stresses on the ice. Winds produce extensive 
deformation, as is evidenced by the occurrence of leads and pressure 
ridges . According to Wittmann and Schule (1966) , up to 10 percent of 
the surface in the Central Arctic consists of open leads during all 
seasons. Untersteiner (1966), however, pointed out that this esti- 
mate is incompatible with several other estimates of the heat and ice 
budget and will require further investigation. Changes in the heat 
budget resulting from the variation of open water areas are not cov- 
ered by the model. 

Wittmann and Schule also report that 13 to 18 percent of the ice 
surveyed was made up of pressure ridges with a thickness of 20 to 30 
meters. Often zones of several hundred square kilometers were more 
than 50 percent covered by "pressure ice." It has been stressed that 
this model computes ice thickness only as a result of thermodynamic 
processes. If mechanical deformation adds a significant contribution 
to the amount of ice stored in the Arctic, then the predictions from 
the model can be used in mass budget studies only in conjunction with 
suitable assumptions concerning the volume of ridged ice. 

An indirect result of the prevailing wind patterns is the removal 
of large quantities of ice by surface currents. As mentioned in the 
introduction, an amount of ice roughly equivalent to a slab 30 cm 
thick and covering an area equal to the whole Arctic Ocean, is ex- 
ported annually. It would be relatively simple to incorporate an 
export factor into the equations. However, the drift patterns indi- 
cate only a slight divergence in the Pacific Gyral, implying that most 
of the exported ice originates in the peripheral seas. Thus the con- 
tribution from the Central Arctic is probably small and has been ne- 
glected. 

Another serious problem occurs in the treatment of the turbulent 
fluxes at the boundaries. In Eqs . (9) and (10), F g , F^, and F w are 
assumed to be independent of growth rate and the physical state of 
the ice or snow at the boundaries. This assumption is not realistic 
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for thin or rapidly growing ice and, ideally, these fluxes should re- 
flect conditions at both surfaces. But unless theoretical models of 
the atmosphere and ocean can be incorporated into the boundary condi- 
tions, the ice must remain isolated from its environment. Although 
the turbulent fluxes are small relative to the other energy fluxes, 
they are not small in comparison with the net energy balance at the 
surface. A series of experiments to be described in Section 4.3 dem- 
onstrates the extent to which these fluxes can influence ice thick- 
ness and temperature. 

Other conceptual omissions include precise treatment of the down- 
ward heat transport from melting snow, the shape of the salinity pro- 
file for various ice conditions, and the heat storage by melt ponds. 
The techniques used to handle the first two problems have been pre- 
viously discussed. No provision has been included to account for the 
heat released by freezing of the melt ponds in the autumn. The com- 
plete freezing of a 40-an-deep melt pond releases nearly 3 kcal/cm^. 
An addition of heat of this magnitude could significantly increase 
upper level temperatures, decrease bottom accretion, and prolong the 
period of bottom ablation; the extent of these effects would depend 
upon the length of time it takes for freezing to occur, along with 
the area covered by ponds and their depth at the onset of freeze-up. 

Finally, it should be stressed that the model cannot predict 
whether the ice would return once it had vanished. In fact, results 
for pack ice thinner than 50 cm would probably be inaccurate, since 
such thin ice is mechanically weak and may show different drift and 
deformation patterns, along with a loss of horizontal uniformity. In 
addition, the assumptions regarding salinity and the turbulent fluxes 
are probably invalid for young, rapidly growing ice. 

Although mechanical motions influence the volume of ice contained 
in the Arctic, they affect ice growth only indirectly. The factors of 
principal concern are those which directly influence the energy budget 
(turbulent fluxes, ponds, etc.). The consequences of these latter 
restrictions are examined in Chapters IV and V. 
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2.5.2 Physical limitations 

The assumption of a constant oceanic heat flux is probably valid 
for the major part of the year; however the stable, low-salinity layer 
formed during the summer could considerably reduce F^. Specifying F w 
as a function of time, while retaining the same yearly total, could 
affect the monthly details of mass changes at the bottom, but it seems 
improbable that a time-dependent F^ would seriously change the mean 
ice thickness . 

Another limitation arises from the coarseness of the input data. 
Energy-flux data are given as monthly averages, but integration of the 
model depends upon specifying these fluxes at daily or hourly intervals. 
The monthly averages must therefore be broken up into smaller incre- 
ments by a smoothing process. Consequently the results will reflect 
the smoothing, and short-term trends or reversals will be absent. Al- 
though the model does predict temperature and thickness on a daily or 
hourly basis, only the monthly values can be accepted with any degree 
of confidence. The same is true with regard to albedo, although here 
the boundary condition has been constructed so that albedo values de- 
pend upon surface conditions to some degree. These limitations are 
not inherent in the model, and more detail can be incorporated when 
the empirical information becomes available. 
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III. METHODS OF SOLUTION 

Solution of the theoretical model discussed in the previous chap- 
ter may be attempted by means of an analog computer. Much of the 
groundwork for such an approach has already been laid (Schwerdtf eger , 
1964; U. S. Naval Civil Engineering Lab ., 1965), but the expense of 
designing and constructing such a system has apparently prohibited 
further work. The alternative is to express the modeling equations 
in finite-difference form for solution on a high-speed digital com- 
puter (IBM 7090/7094). The latter method has been chosen for the 
present study. 

3.1 DIFFUSION EQUATION 

Because the heat-conduction equation is parabolic, its solution 
involves specifying a temperature profile at some initial time, t = 0, 
and values at both boundaries during each succeeding time step. The 
problem is then to predict future temperature profiles on the basis of 
this information. Initial -value problems of this sort are tradition- 
ally solved by either an explicit or implicit method. Unfortunately, 
implicit methods are inapplicable here because the equations require 
moving boundaries . 

In choosing a method of solution, the prime consideration must 
be to minimize computer time , while making accurate approximations . 
Untersteiner (1964) solved Eq. (5) by employing the standard forward- 
difference method (Richtmyer , 1957; Forsythe , 1960), which required a 
time step (At) of one hour for a corresponding depth increment (Az) 
of 10 cm; 45 minutes of computer (IBM 709) time were required for only 
one year of integration. The present model Is much more complex and 
it is anticipated that 30 to 50 years of integration will be necessary 
for the ice to achieve equilibrium. These facts, combined with the 
need to integrate at least 25 separate cases, imply that the use of 
the forward— difference method would require an excessive amount of 
time on this computer. 

Since 10 cm was deemed to be the minimum depth resolution nec- 
essary for meaningful results, various explicit schemes were studied 
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and tried, but with little success. Finally, a technique suggested 
by Sauliev (1957a, 1957b) was found to be a powerful method which com- 
bines some of the advantages of both implicit and explicit schemes, 
but is in fact explicit. Because it is unconditionally stable, the 
Sauliev method does not restrict the choice of grid spacing and is, 
in addition, ideally suited to handling nonlinear terms and moving 
boundaries . 

The Sauliev method differs from other explicit methods by ex- 
pressing the second-order differential in terms of values in two ad- 
jacent time levels. The scheme thus appears to be implicit; however 
the implicit terms are asymmetrically arranged so that, if solution 
is initiated at a known boundary value, each unknown can be obtained 
sequentially. Although the Sauliev method was originally designed to 
solve the one-dimensional diffusion equation, Larkin (1964) has suc- 
cessfully extended it to a two-dimensional case and it has been ap- 
plied to other equations as well ( Holton , 1967). 

For the present case, let z and t space be subdivided into J and 
M sub intervals of lengths Az and At, respectively, containing points 
(z,t). Let us now consider Eq. (6) (neglecting the last term for rea- 
sons which will become obvious below). Applying the Sauliev method, 
Eq. (6) has the finite-difference form (see Fig. 4a): 



,t) 



k i [ *T(z+Az.t)-T(z,t) _ f(z.t+At)-T( z-Az,t+At) 



-£ — [T(z+Az,t) + T(z-Az,t+At)] 



where 

k. At 
x 
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Equation (16) may be solved only if T(z-Az , t+At) is known. If z = 2, 
then T(z-Az,t+At) is a boundary value and T(2,t+At) may be evaluated, 
which in turn allows solution of T (3, t+At) and so forth, as z ascends 
from 2 to J-l. The effect of a boundary value is thus propagated 
"through the entire matrix at a given time level, which is not the 
' case in other explicit methods. It may be noticed that all values of 

T(z,t+At) may be determined from only one set of boundary values; 

however, equivalent values could have been determined from the other 

boundary (z - J) . The formulation of this alternative estimate is 

completely analogous to (16) : 



where T refers to values obtained as z descends from J-l to 2 (see 
Fig. 4b). The Sauliev method is unique in that it offers two indepen- 
dent estimates of the temperature at (z,t+At). This fact offers an 
opportunity to reduce one of the more serious limitations of the meth- 
od. 

Although it may easily be shown that (16) and (17) are stable re- 
gardless of the choice of At and Az, the truncation error associated 
with the Sauliev method is an order of magnitude greater than with the 
forward-difference scheme. Taylor Series expansions may be used to 
determine truncation errors in (16) and (17) (Larkin, 1964). Examina- 
tion of these expansions show that alternate terms in the two equa- 
tions are equal and opposite, which immediately suggests the possibil- 
ity of combining the expansions to reduce the truncation error. If 
(16) and (17) are averaged, we have: 

5(Zft+At) - = T(z,t+At) +T(z,t+At) _ = ^1 T(Zf t) + 

- T(z-Az,t) + T(z+Az,t) + T(z-Az,t+At)] . (18) 



Calculation of T(z,t+At) 



T(z+Az,t+At) 
b. Calculation of T(z,t+At) 



Fig. 4 -- Calculations by the Sauliev Method <, 
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The result of averaging T and T is that two of the three largest error 
terms drop out, leaving a truncation error for f roughly proportional 
to At 2 ; this is comparable to other methods. The increase in accuracy 
justifies the added computer time necessary to evaluate T(z,t+At). 

It is of interest to note that Eq. (18) has the exact form of the 
well known Crank-Nicholson method (Dingle, 1965). However, T(z+Az,t+At) 
and T(z-Az,t+At) do not result from simultaneous solutions; rather 
they are approximated by the previously discussed methods. Equation 
(18) is thus an explicit approximation to an implicit scheme. 

In dealing with the penetration of short-wave radiation, it is 
necessary to treat the I q term in Eq. (6) as a forcing function be- 
cause of the unique form of the Sauliev equations . Assume that z = 1 
is the upper surface, and consider the situation at time t if the I q 
term were included in (17). Absorption of solar energy would cause 
T(2,t+At) to have a value higher than from pure conduction, but 
T(2',t+At) is used to determine T(3,t+At). T(3,t+At) thus would feel 
the'effect of solar heating at T(2,t+At), in addition to direct pene- 
tration of short-wave radiation to this level. Since no conduction 
may occur during a particular instant of time, the described situa- 
tion is opposed to reality and hence f must first be evaluated and 
1 the I q term treated as a forcing function: 

f (Zf t+At) = IkJ **t) + T(z,t + At) + liVl exp ( _ <iZ) . (19 ) 

Equations (3), (4), (16), (17), and (19) are used to calculate temper- 
atures within the ice. If snow is present, I q = 0. 

As a test of the Sauliev method, Unters teiner ' s (1964) experiment 
was repeated. Equation (6), the simplified form of Eq. (5), was inte- 
grated using the boundary values and grid spacing of Untersteiner . 
The maximum deviation of the results from those obtained by the for- 
ward-difference method was less than 0.05°K. This demonstrates the 
validity, not only of the technique, but also of the assumptions made 
to obtain (6) . 
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The stability condition for the forward-difference method, 
E, < 1/2, approximates the criteria for most of the other explicit 



method was stable, but for large values of ? max the truncation errors 
became large and the resulting temperature field lost most of its de- 
tail. On the basis of these experiments, it appears that truncation 
errors from the Sauliev method are unimportant if C max is less than 
8. For the initial test of the model, At = 12 hours and Az = 
10 cm (5 * 6) were chosen. With this choice of grid spacing, a 
yearly temperature field may be calculated from (6) , with about 15 
seconds of IBM 7094 computer time and a deviation from the forward- 
difference method of less than 0.1°K(0.05 percent). 



3.2 BOUNDARY EQUATIONS 

The balance-of-energy equations , together with the heat-conduc- 
tion equation, give surface temperatures and mass changes at the bound- 
aries. The greatest difficulty in a numerical solution to these equa- 
tions is the moving boundaries. The coordinate system is fixed in 
space, but has two floating grid points located at the boundaries. 
As the ice grows and shrinks, the slab moves up and down through the 
grid system. Time and depth increments are constant, except for the 
space interval next to either boundary, which is some fraction of Az. 
The only assumption necessary for this treatment is that the tempera- 
ture profile is linear within 2Az of the surfaces. This assumption 
is valid at the bottom, where the temperature gradient is nearly lin- 
ear throughout the year, and likewise at the upper boundary during the 
ablation season. However the linear assumption must also be made dur- 
ing the first 2Az of snow accumulation and this could lead in the model 
to heat losses slightly greater than reality. 

3.2.1 Ice/water interface 

The simplest situation exists at the bottom of the ice. Assume 
that, with some initial temperature profile at time t, the bottom of 
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the ice lies a distance Az' below the nearest grid point. At time 
t + At there will exist a new profile and the bottom will have moved 
a distance Ae (see Fig. 5a). If a linear temperature is now assumed 
to exist between T(J,t + At) and T(J - 2,t + At), then: 

T(J-l,t+At) = T(J,t+At) + A^t'+L t^ J - 2 » t+&t > - T(J,t+At>] . (20) 



Applying Eq. (17) at (J-2,t+At): 

(l+£)T(J-2,t+At) = (l-C)T(J-2,t) + £[T(J-l,t+At) +T(J-3,t)] . (21) 

Substituting (20) into (21) : 

£AzT(J. t+At)+(A2+Az'+Ae)r(l-C)T(J-2,t)+CT(J-3,t)] 
T(J-2,t+At) = (1+OAz + Az' + Ae 

(22) 

Equation (6), the lower-boundary balance equation, has the finite- 
difference form: 



k.[T(J,t+At) - f(J-2,t+At)] 

F — — 

w Az + Az' + Ae 



Incorporating (22) into (23) yields: 

Ac _ l| F w , (1+QAz+Az'l l fd+QAz+Az 
At " 2lq b + At J + At 

k. 1 f 

— i- [T(J,t+At) + (5-DT(J-2,t) - ?T(J-3,t)] . (24) 
q b At j 

The right-hand side of Eq. (24) is written in terms of known quanti- 
ties and we may now solve for Ae . Knowing Ae , it is then possible to 
solve for T(J-2,t+At) from (22) and for T(J-l,t+At) from (20). 




T(J-3,t) 





T(4,t) 1 . 

h At »-| 

At the upper boundary during surface ablation 
Fig. 5 -- Geometry of the grid system. 
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3.2.2 Upper surface 

Because of the complexity of the balance equation and the many 
processes occurring, the solution of temperature and mass changes at 
the upper boundary offers more difficulties than at the lower. Assume 
that the upper surface is located at z = 1 and that no ablation or 
accretion is taking place, i.e., T(l,t+At) < 272. 9°K. According to 
(16): 

(l+O T (2,t+At) = (l-O T (2,t) + C[T(l,t+At) + T(3,t)] . (25) 
Combining (25) with the finite-difference form of Eq. (8): 



a k MC 1 -^) T (2»t) + £T(3,t)] 
e T aT(l,t+At) H + — T(l,t+At) - — 



Iq _ (1 _ a) F r - F L - F a - F £ - 0 . (26; 

Equation (26) is a 4th order equation in T(l,t+At) and is solved by 
the Newton-Raphson iterative method, a recurrence formula (Hildebrand, 
1956). 

If the surface temperature, T(l,t+At), is below the freezing 
point, Eq. (25) is solved and Eq. (19) is then used to obtain 
T(z,t+At). If, on the other hand, (26) predicts that the surface tem- 
perature is above the freezing point, then melting must occur and the 
boundary will move. Assuming a linear temperature profile within 2Az 
of the surface and using the same notation as we did for the bottom 
of the ice (see Fig. 5b), we have: 



Az'+ Ae 

T(2,t+At) = T(l,t+At) + — [T(3,t+At) - T(l,t+At)] . (27) 

Az+Az'+Ae 
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Applying (16) at (3,t+At) and including (27): 



£AzT(l,t+At) + (Az+Az'+Ae ) [ (l-£)T(3,t)+£T(4,t)] 

T(3,t+At) = - . (28) 

(l+C)Az = Az' + Ae 



If Ae Q ^ 0, Eqs. (8) and (28) give: 
K 



e oT(l,t+At) + - 



(l+£) Az+Az'+Ae 



[T(l,t+At) + (C-l)T(3,t) - £T(4,t)] 



I - (1-a) F - F - F - F = q 
o r L s I 



Solving (29) for Ae : 

Ae o e L aT(l,t+At) 4 + I Q - (l-a)F^-F T -F c> -F o 
A t 2 % 

L (l+5)Az+Az' £ aT(l,t+At) 4 + I - (l-a)F -F -F -F 



[T(l,t+At) + (?-l)T(3,t) - ^T(4,t)]> 



Since Ae * 0, T(l,t+At) is known to be at the freezing point of the 
ice or snow. If Ae^ is positive, then ablation has ceased and 
T(l,t+At) is again determined from (26). After the snow has become 
isothermal, Eq. (30) reduces simply to: 

^ = [e L oT(l,t+At) 4 + I Q - (l-a)F r - * L - *, - \]/ % . (3: 



Simultaneous solution of Eqs. (3), (4), (13), (16), (17), (20), (22), 
(24) , (27) , (28) , and (30) or (26) allows the complete determination 
of T(z,t+At), f(z,t+At) and any mass changes that may be taking place. 
Finally, f(z,t+At) is calculated from Eq. (19). 

Actually, solution of the model is not quite this simple during 
the snow-accumulation period. If the snow depth is in excess of 2Az, 
then Eqs. (27), (28), and (29) may be applied, specifying that Ac q = 0. 
If the snow accumulation is less than 2Az , a linear temperature profile 
in the entire snow layer must be assumed and the ice/snow boundary 
condition incorporated into equations analogous to (27) , (28) , and 
(29) . Because a second order differential cannot be accurately evalu- 
ated across dissimilar media, special techniques must be used to find 
the temperature at the ice/snow interface.^" Using these methods, cal- 
culation of an annual temperature field and the associated thickness 
changes for 3-meter ice requires only 20 sec of IBM 7094 computer time 
per year of real time. Subroutine YARIT, a FORTRAN IV description of 
this integration, is presented in the appendix. 

3.3 SMOOTHING THE INPUT ENERGY FLUXES 

Since values of the energy fluxes are available only in the form 
of monthly averages, some method is needed to realistically divide 
them into increments corresponding to the chosen time steps. An ana- 
lytic expression describing energy variations with time is the most 
convenient method. To attain smoothed values, it was assumed that the 
monthly averages were representative of the 15th day of the month. 
The monthly values were then divided by the number of time increments 
in a month and a curve of energy vs. time constructed. The irregular- 
ity of the curve precluded any attempt to describe the entire curve; 
however portions of it were approximated by polynomials. Since any 

"'"Briefly, T(z,t+At) is found in the snow and T(z,t+At) in the ice. 

Temperature at the interface is calculated from (12) , which in turn 

allows us to determine T(z,t+At) in the ice and T(z,t+At) in the snow. 
T(z,t+At) is then evaluated from (19). 
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n + 1 points determine an nth-order polynomial, four points on the 
curve were used to find a 3rd-order polynomial of the form 

F(t) = A x t 3 + A 2 t 2 + A 3 t + A 4 . (32) 

The term F(t) directly gives the energy flux in a time step At, cen- 
tered around time t. Values for F(t) are calculated only between the 
two central points on the energy vs. time curve, and smooth trends can 
be established from month to month. The constants in (32) are reeval- 
uated each month to produce 12 different polynomials. This method 
has the disadvantage that the monthly sums of the flux increments may 
deviate from the specified amount by as much as 0.5 percent. Subrou- 
tine FLIP accomplishes this smoothing. 

3.4 ANALYTIC EXPRESSION FOR ICE SALINITY 

As discussed in Section 2.4.3, Schwarzacher ' s salinity profile 
may be approximated by an equation of the form: 



S(z,t) = A + : 



The constants in (13) may be evaluated by imposing the following con- 
ditions: the salinity of ice at the boundaries z = 0 and z = H is 
fixed and has values S(0,t) = S q and S(H,t) = S R ; the slope of the 
profile is zero at both boundaries; and finally the salinity has val- 
ues and S 2 at z = z^' and z = z^ x . Both z^' and are some con- 
stant fraction of H; i.e., z^ = r^H and z^ = n 2 H, where n ] _ and n 2 
are less than unity. Applying these constraints and solving for the 
constants, Eq. (13) becomes: 



\ • 
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2S,-S n -SJ 



How well (33) fits the observed profile depends upon the choice of 
n 1 and n 2 - In the present case, - 0.3 and n 2 = 0.6 produced the 
best fit. Using these values, the maximum deviation of (13) from 
Schwarzacher's profile was less than 5 percent. Subroutine SALPR 
calculates the salinity profile in the program. 
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IV. RESULTS OF SELECTED INTEGRATIONS 

Arbitrary initial conditions of temperature and ice thickness are 
imposed upon the equations, and integration proceeds until equilibrium 
is reached, that is, by our definition, when the annual bottom accre- 
tion is within 1 millimeter of the surface ice ablation. Steady-state 
solutions are generally reached after 30 to 100 years. Although there 
is a continuing output of data during the integration, usually only 
the steady-state annual patterns of temperature and thickness will be 
shown in the text. The isotherms are calculated by a contour plotting 
subroutine, and the temperature field is graphed by a Calcomp plotter. 
Each graph is accompanied by a table containing monthly averages and 
totals of various pertinent results. 

4.1 INITIAL TESTS OF THE MODEL 

To check the validity of the model and the numerical approxima- 
tions employed, the best available estimates of the present energy 
fluxes were applied to the model. Assuming that arctic pack Ice is 
in near-equilibrium, predictions of ice thickness and temperatures 
from this test were compared with Soviet and American observations. 
Investigations of response time and salinity effects were then con- 
ducted. 

4.1.1 Heat budget of Fletcher 

With the heat budget suggested by Fletcher (Table 1), with the 
other input data discussed in Chapter II, and with an initial thick- 
ness of 340 cm, equilibrium was reached after 38 years. Figure 6 
illustrates the steady-state temperature and thickness patterns for 
this case. In Fig. 6, and all the similar plots that follow, time is 
graphed along the abscissa and thickness along the ordinate. The me- 
dia (air, ice, and water) are distinguished by shading. Boundaries 
of the ice/snow slab appear as heavy black lines and isopleths of tem- 
perature as thin solid curves. The isotherms in the ice are labeled 
in negative degrees Celsius; Isotherms in the snow (not labeled) are 
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drawn at 2°C intervals. Because hydrostatic adjustment of the ice has 
been removed from the equations , the upper and lower boundaries move 
independently of one another. Thus, instead of the upper boundary oc- 
cupying a nearly constant level and the lower surface showing large 
changes in depth, the slab migrates downward with time. The advantage 
is that mass changes at each boundary are separately indicated, but 
it must be remembered that the depth indicated on the ordinate does 
not necessarily refer to ice thickness. Table 2 gives monthly and 
yearly values of ice thickness, surface temperatures, net radiative 
fluxes, and the conductive fluxes at both boundaries; also included 
in the table are values of the annual mass changes and associated 
dates . 

From Table 2 and Fig. 6, it may be seen that the model predicts 
an average equilibrium thickness of 288 cm, with a maximum of 314 cm 
and a minimum of 271 cm. It is difficult to assess how representative 
these values are. Early submarine data ( Lyon , 1961) indicate an aver- 
age thickness of 4 to 5 meters, but this includes pressure ridges and 
may be biased by the location of the submarine's track. Later sub- 
marine observations ( Wittmann and Schule , 1966) show 2 meters to be 
the most frequent ice thickness, but again, these data are regionally 
limited. Most of the drifting stations were established on pack ice, 
which averaged close to 3 meters in thickness ( Petrov , 1954). Typical 
results are those of Untersteiner (1961) who found that thickness var- 
ied between 315 and 250 cm during the year. Considering that the in- 
put data are derived from these drifting stations, it is to be expec- 
ted that the predicted thicknesses should agree most closely with 
their observations. 

Perhaps more indicative is the pattern of mass changes predicted 
by the model. According to Soviet data ( Yanes , 1966), average abla- 
tion on polar ice is 37 cm; snow melt begins in the first half of June 
and ablation generally ends between 10 and 23 August. An average from 
all U. S. drifting stations from 1957 to 1963 gives a surface ablation 
of 42 cm (Hanson , 1965). Table 2 shows a predicted ice ablation of 
40 cm, starting on 29 June and terminating on 19 August. Snow melt 
begins on 8 June. All these features correspond closely with obser- 
vations. 
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At the lower surface there are 45 cm of accretion and 5 cm of abla- 
tion. Unfortunately, mass changes at the bottom have never been well 
determined in the field. Yanes (op. cit.) feels that bottom ablation 
in equilibrium ice is small, if it occurs at all. Hanson (1965) mea- 
sured a value of 10 cm, and Untersteiner (1961) observed 20 cm of ab- 
lation and 50 cm of accretion. In view of these uncertainties, the 
predicted values are not unreasonable, but lack of empirical data makes 
it impossible to evaluate the timing of these changes accurately. 

Temperature profiles provide another basis for evaluating the 
theoretical results. An observed temperature field ( Untersteiner , 
1961) is reproduced in Fig. 7. Comparison of Figs. 6 and 7 shows good 
agreement throughout, except for fall temperatures within the ice. 
Smoothing of the input data is reflected in the results of Fig. 6 and 
the usual January warming trend ( Laktionov , 1955; Hisdal, 1960) does 
not appear. Figure 7 suggests that ice temperatures in the theoreti- 
cal model respond to surface conditions more rapidly than do the tem- 
peratures in real ice. Neglecting latent-heat release from surface 
melt ponds, underestimating I , specifying linear heat loss through 
the snow layer, or improper assumptions regarding the salinity pro- 
file could all contribute to this discrepancy. 

Table 3 shows monthly surface temperatures from numerous drifting 
stations, averaged over 26 years of observations ( A. Hanson , personal 
communication), and a comparison with the theoretical results. The 
annual averages are remarkably close, but temperatures predicted by 
the model are slightly colder in the fall and warmer in the spring 
than the observed ones. It should be remembered that the calculated 
temperatures apply for the surface and a comparison with observed air 
(screen) temperatures would have to take into account that, generally, 
a ground inversion prevails from September to April, whereas in May 
lapse conditions prevail. We suspect that the differences shown in 
Table 3 are more or less accidental, as some are positive and some 
are negative. 

The net radiative fluxes calculated by the model also agree fa- 
vorably with observations. According to Marshunova (1961) and Fletcher 
(1965), the net short-wave radiation ranges from 18.6 kcal/cm^ year 
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Table 2. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS 



Symbol 


Variab le 




Feb 


Mar 




T s o 


Mean snow surface 
temperature ( C) 


-30.98 


-33.12 


-31.85 


-22.46 


T i o 


Mean ice surface 
temperature ( C) 


-17.95 


-19.04 


-19.18 


-15.84 


H 


Mean ice thickness (cm) 


281.7 


289.0 


296.6 


304.1 


F 

c ,o 


Heat flux through surface 
(kcal/cm 2 ) 


0.781 


0.823 


0.716 


0.355 


F c,h+H 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.580 


-0.605 


-0.623 


-0.574 


(l-Qf>F r 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


4 

F L _aT o 


Net long-wave radiation 
(kcal/cm 2 ) 


-1.967 


-1.562 


-1.814 


-2.451 


CM 




JG | SEP | OCT 


,..J..* 3v ...J.. 0f - 








Fig. 6 -- Annual equilibrium temperature and thickness field 
for Case 1, Fletcher's heat budget. 



IN CASE 1, FLETCHER'S HEAT BUDGET (SEE TABLE 1). 



May Jun Jul Aug Sep Oct Nov Dec Year 



-9.44 


-0.33 


-0.10 


-1.12 


-10.04 


-20.37 


-28.74 


-30.28 


-18.24 


-9.71 


-2.31 


-0.10 


-0.85 


-6.24 


-10.70 


-14.84 


-17.16 


-11.16 


310.3 


314.0 


296.1 


278.1 


272.8 


270.8 


270.7 


275.0 


288.3 


-0.022 


-0.119 


-0.042 


0.247 


0.663 


0.781 


0.884 


0.807 


5.874 


-0.448 


-0.290 


-0.135 


-0.036 


0.004 


-0.015 


-0.257 


-0.510 


-4.069 


3.129 


4.395 


4.959 


2.710 


0.605 


0.085 


0 


0 


18.181 


-2.200 


-1.935 


-0.903 


-0.984 


-0.712 


-0.781 


-1.438 


-1.592 


-18.339 



Date Onset End 

8 June Snow ablation 

29 June Ice ablation 

17 July Bottom accretion 

19 August Surface ice ablatio 

4 November Bottom accretion 



Other Data 

Surface ice ablation* 40.1 cm 
Bottom ablation 5.1 cm 

Bottom accretion 45.2 cm 
Short-wave radiation penetrating 

surface during snow-free period 1.291 kcal/cm 

Equals net bottom accretion. 




Fig. 7 -- Observed temperature of perennial ice at IGY Station Alpha, 
1957-1958 (after Untersteiner , 1961). 
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at 75°N to 16.9 at 90°N, while net long-wave radiation ranges from 
-21 to -19.4 kcal/cm^ year. Solar radiation absorbed by the surface 
in the model agrees closely with empirical data, implying that the 
albedo assumptions were adequate; however, the long-wave loss is small- 
er than observations indicate, primarily during the winter months. 
This smaller loss appears to arise from two causes. First, the model 
deals only with thick ice covered with a uniform blanket of snow; ob- 
servations in the field are made over ice with variable thickness and 
snow cover. Second, necessary smoothing of the input functions re- 
moves short-term fluctuations in the surface temperature. Because of 
the nonlinear form of the Stef an-Boltzmann Law, this produces a bias 
toward lower surface temperatures and hence a smaller long-wave loss. 

In a state of equilibrium, the sum of all energy fluxes into the 
ice should balance the sum of those leaving it. The results in Table 
2 show that about 0.5 kcal/cm^ year more energy leaves the slab than 
enters it. However, there are three energy sources which are not seen 
explicitly in the output data: 

(a) Short-wave radiation that penetrates the ice and is not 
absorbed, passing instead through the ice and into the ocean. This 
flux is small for thick ice (25 cal/cm year in the present case) , 
but for thin ice it may reach values in excess of 25 percent of I . 

(b) Latent heat released within the ice. As sea water freezes 
at the bottom of the ice, small pockets of brine are trapped. Con- 
tinued accretion moves the older ice upward into regions of lower 
temperature, causing freezing in the brine pockets and a reduction 
in brine volume. The resulting heat release is accounted for by 
changes in (pc^ [see Eq. (3)]. 

(c) When the snow is saturated with melt water and is isother- 
mal, no heat flows toward the surface. The ice/snow interface, how- 
ever, remains at the freezing point, causing a large conductive flux 
into the ice. The sudden introduction of this flux is responsible 
for the step seen in the -1.0, -0.5, and -0.25 isotherms during June 
(Fig. 6). 
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When the effects of these energy sources are included, the energy bal- 
ance for the entire slab presumably becomes zero and there is no con- 
tradiction. 

Despite the many uncertainties in the input data, the picture 
presented by the model is surprisingly consistent with field observa- 
tions. Thus the model and the various approximations appear to be 
valid and lend confidence to the model's application to other possible 
situations. To facilitate reference, all cases computed, and their 
characteristic difference for input data from Case 1, are listed in 
Table 29 on p. 130. 

4.1.2 Response time 

A study of the nonequilibrium behavior of sea ice offers insight 
into the nature of ice growth. To determine the response of ice 
thickness to a change of the external forcing, we assume that the ice 
is initially at equilibrium under the "standard" conditions of Case 
1. A step change in F w from 1.5 to 4.5 kcal/cm 2 year is then imposed. 
The thickness on 1 January for each succeeding year is graphed in 
Fig. 8. The exponential shape of this curve is typical of the manner 
in which both temperature and thickness approach equilibrium and is 
independent of the specified boundary conditions. The curve resembles 
that for the discharge rate of a simple capacitor. This similarity 
suggests that temperature and thickness response might be simulated 
by a relatively simple electrical analogy; for a given capacitor, the 
size of the initial charge would correspond to the deviation of tem- 
perature or thickness from equilibrium values. If this capacitor anal- 
ogy is valid, the response time (period from change of input to achieve- 
ment of equilibrium) of sea ice would be independent of the initial 
thickness. To check this, the experiment was repeated starting from 
initial thicknesses that differed by 10 to 400 cm from the equilibrium 
value. Results are shown in Fig. 9. There is little variation in re- 
sponse time between large deviations. However, when the initial thick- 
nesses are close to the equilibrium value, there is a rapid variation 
in response time and the capacitor analogy breaks down. Although a 
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comprehensive series of tests was not run for cases with greater equi- 
librium thicknesses, various isolated cases indicate that, as equilib- 
rium thickness increases, the region of rapid variation in response 
time becomes smaller. 

Despite the limited dependence on initial thickness, response 
time is determined primarily by the equilibrium thickness . As thick- 
ness increases, the ice responds more and more sluggishly to surface 
conditions. Mass and temperature changes at the upper boundary are 
not greatly influenced by temperature gradients within the ice and 
snow, but mass changes at the bottom are directly determined by the 
magnitude of these gradients. Therefore, even though surface abla- 
tion rapidly achieves a constant value, net bottom accretion continues 
to change slowly. It is clear from these considerations that the re- 
sponse time must increase with equilibrium thickness. 

These results show that the response time is asymmetric. For 
example, when is increased from 1.5 to 4.5 kcal/cm 2 year, the re- 
sponse time is about 20 years; if the ice is at equilibrium for an F 

2 2 w 

of 4.5 kcal/cm year and F w is decreased to 1.5 kcal/cm year, the re- 
sponse time is 44 years. An anomalous year will produce nonequilib- 
rium conditions of several years duration, because of the exponential 
shape of the growth curve. Because response time is principally de- 
termined by the final thickness, the magnitude of the anomaly plays 
only a small part in determining the duration of the nonequilibrium 
conditions . 

4.1.3 Effects of ice salinity 

The salinity profile used in the model (Fig. 3) is probably char- 
acteristic of summer conditions only. However, it was assumed that 
any deviations from the profile were relatively unimportant during 
the winter when the thermal constants of the ice are nearly indepen- 
dent of salinity. In an attempt to verify this assumption and to de- 
fine more precisely the effects of ice salinity, we have assumed in 
Case 2 that the salinity has a value of 0.09°/ o6 and is constant 
throughout the slab. Figure 10 shows the steady-state temperature 
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Fig. 9 - 



Response time of sea ice as a function of the deviation of 
the initial thickness from equilibrium thickness, F w = 
4.5 kcal/cm 2 year. 



Table 4. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 2, 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 


-31.17 


-33.30 


-32.02 


-22.60 


-9.51 


Mean ice surface 
temperature ( C) 


-18.80 


-19.82 


-19.95 


-16.52 


-10.10 


Mean ice thickness (cm) 


306.1 


313.7 


321.5 


329.1 


334.9 


Heat flux through surface 
(kcal/cm 2 ) 


0.742 


0.788 


0.682 


0.325 


-0.039 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.666 


-0.680 


-0.692 


-0.629 


-0.454 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-1.927 


-1.526 


-1.780 


-2.420 


-2.183 



JUL , AUG | SEP | OCT | NOV | DEC 




Fig. 10 - 



•- Annual equilibrium temperature and thickness field 
for Case 2, 2 constant salinity (0.09 °/ 00 ), F w = 
1.5 kcal/cm year. 



CONSTANT SALINITY (0.09 °/oo). ALL OTHER INPUT DATA AS IN CASE 1. 



Jun Jul Aug Sep Oct Nov Dec Year 



-0.34 


-0.10 


-0.93 


-9.62 


-20.62 


-29.06 


-30.53 


-18.32 


-2.32 


-0.10 


-0.59 


-5.35 


-11.73 


-16.21 


-18.22 


-11.64 


337.8 


317.3 


294.7 


289.1 


288.4 


292.3 


298.8 


310.3 


-0.115 


0.011 


0.331. 


0.775 


0.722 


0.818 


0.756 


5.795 


-0.212 


0.051 


0.127 


0.082 


-0.265 


-0.515 


-0.639 


-4.491 


4.395 


4.959 


2.710 


0.605 


0.085 


0 


0 


18.181 


-1.934 


-0.903 


-1.037 


-0.825 


-0.722 


-1.372 


-1.541 


-18.170 



Date 

8 June 
24 June 
28 June 
20 August 

4 October 



Onset 
Snow ablation 

Ice ablation 

Bottom accretion 



Bottom accretion 



Surface ice ablation 



Other Data 

Surface ice ablation* 

Bottom ablation 

Bottom accretion 

Short-wave radiation penetrating 
surface during snow- free period 



41.5 cm 
8.8 cm 
50.3 cm 

1.291 kcal/cm 2 



Equals net bottom accretion. 



-64- 



and thickness patterns. Comparing these results with those from Case 
1, we see that ice thickness is greater by 20 to 25 cm the year around, 
and that ice temperatures are 1 to 2°C lower during fall and winter. 
The greatest difference occurs, as expected, during the summer. Be- 
cause of the negligible brine volume, summer heating (see -1°C isotherm) 
penetrates to much greater depths. The steepness of the -2°C isotherm 
shows that the entire slab responds rapidly to temperature changes at 
the surface. 

The temperature pattern in Case 2 shows the extent to which the 
brine pockets store heat in the summer and retard heat transport be- 
tween the boundaries. Reduction of salinity near the bottom of the 
ice allows heat to be more efficiently conducted upward into the ice 
and thus increases bottom accretion. At the same time, the larger con- 
ductive heat fluxes also cause slightly greater ice ablation at both 
boundaries. In this case, the increased bottom accretion is the dom- 
inant effect, producing the predicted thickening. 

Case 3 is similar to Case 2, except that F w was increased to 
4.5 kcal/cm year. Figure 11 and Table 5 give the results of the in- 
tegration. Comparison of Cases 3 and 1 confirms the conclusions reach- 
ed in Case 2 and also shows that the thickness of the ice does not ap- 
preciably change the effects of salinity. 

The model may also be used to compare the growth of ice in fresh 
water and sea water. In Case 4, we have again assumed that the salin- 
ity is a constant 0.09°/ oo and that the ice temperature at the lower 
boundary is -0.1°C; i.e., the slab is floating in essentially fresh 
water. All independent energy fluxes are the same as in Case 1. Re- 
sults are presented in Fig. 12 and Table 6. The average equilibrium 
thickness is 348 cm, 60 cm greater than Case 1 and 40 cm greater than 
Case 2. The major difference between Cases 2 and 4 is the 50 percent 
less bottom ablation of Case 4. This is the result of the steeper 
temperature gradients near the bottom from September through March in 
Case 2 and is responsible for the predicted thickness increase. Al- 
though the assumption regarding heat flux in the water is certainly 
not characteristic of a fresh water body, Case 4 points out the main 
difference between the growth of sea ice and perennial lake ice; fresh 
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ice cannot store large amounts of heat and instead transfers it rapid- 
ly between the boundaries. Consequently, under similar conditions, 
fresh ice grows more rapidly and to a greater equilibrium thickness 
than does sea ice. 

The results of Cases 2 to 4 indicate that the ice salinity has 
little influence on ice thickness but markedly influences summer tem- 
perature profiles. Therefore, since the winter temperature differ- 
ences between Cases 1 and 2 appear to result solely from the slight 
thickness change, it may be assumed that reasonable ice temperatures 
can be obtained from the model knowing only the summer salinity pro- 
file. 

4.1.4 Heat budget of Vowinckel and Orvig 

An obvious application of the sea-ice model lies in its ability 
to evaluate the internal consistency of a given energy budget. We 
have seen in Section 4.1.1 that Fletcher's (1965) heat budget produces 
results in good agreement with observations; however, it seems likely 
that other combinations of input functions could produce equally ac- 
ceptable results. A heat budget markedly different from that of 
Fletcher was presented by Vowinckel and Orvig (1966, 1967). Their 
values for F^, F r , F g , F £ , net short-wave radiation, and implied al- 
bedos are reproduced in Table 7. Their surface energy balance was 
determined according to the following equation (using our notation) : 

(1 - a)F r + F L - e^T* + F g + F £ + 0 = 0 , (34) 

where 0 is the "heat exchange with deeper layers of the ocean." Ob- 
viously 0 is a residual term; it appears to include the effects of 
oceanic heat flux, penetration of short-wave radiation, heat conduc- 
tion to the surface, surface ablation, and mass changes at the bottom 
of the ice. According to Vowinckel and Orvig, the annual sum of. 0 

(5.93 kcal/cm year) corresponds to the oceanic heat flux (F ). if 

w ' 

the ice is in equilibrium. This is not strictly true because the en- 
ergy lost at the surface through snow melt is not regained at the 



Table 5. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 3, CONSTANT SALINITY 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature ( C) 


-29.31 


-31.46 


-30.26 


-21.18 


-8.66 


Mean ice surface 
temperature ( C) 


-10.56 


-11.65 


-11.90 


-9.38 


-4.97 


Mean ice thickness (cm) 


100.0 


110.0 


119.9 


127.7 


130.3 


Heat flux through surface 
(kcal/cm 2 ) 


1.123 


1.156 


1.040 


0.644 


0.183 


Heat flux through bottom 
(kcal/cm 2 ) 


-1.085 


-1.106 


-1.054 


-0.777 


-0.350 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-2.309 


-1.895 


-2.138 


-2.740 


-2.406 


CM 
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Fig. 11 -- Annual equilibrium temperature and thickness field 
for Case 3,-constant salinity (0.09 °/ 00 ), F w = 
4.5 kcal/cm year. 



(0.09 %o), F w " 4.5 kcal/cm 2 year. ALJ, OTHER INPUT DATA AS IN CASE 1. 



Jun 


Jul 


Aug 


Sep 


Oct 


Nov 


Dec 


Year 


-0.23 


-0.10 


-1.36 


-9.75 


-19.47 


-27.40 


-28.63 


--17.32 


-0.84 


-0.10 


-1.16 


-5.39 


-7.19 


-8.94 


-9.92 


-6.83 


126.7 


96.4 


70.3 


63.7 


70.2 


79.3 


89.9 


98.7 


-0.046 


0 


0.116 


0.737 


0.990 


1.169 


1.149 


8.262 


0.064 


0.417 


0.252 


-0.688 


-0.949 


-1.112 


-1.125 


-7.512 


4.507 


4.959 


2.710 


0.605 


0.085 


0 


0 


18.293 


-1.966 


-0.903 


-0.916 


-0.787 


-0.989 


-1.722 


-1.934 


-20.705 



Date Onset End 

13 May Bottom accretion 

6 June Snow ablation 

27 June Ice ablation 

17 August Surface ice ablation 

27 August Bottom accretion 



Other Data 

* 

Surface ice ablation 41.8 cm 

Bottom ablation 26.7 cm 

Bottom accretion 68.5 cm 

Short-wave radiation penetrating 2 
surface during snow-free period 1.326 kcal/cm 



Equals net bottom accretion. 



•68- 



Table 6. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 4, CONSTANT SALINITY 



Variable 


Jan 


Feb 


Mar 


Apr 




Mean snow surf age 


-31.16 


-33.31 








Mean ice surface 


-18.75 


-19 .87 








Mean ice thickness (cm) 


343.3 


350.4 


357.8 


365.1 


371.1 


Heat flux through surface 
(kcal/cm 2 ) 


0.744 


0.785 


0.678 


0.321 


-0.041 


Heat flux through bottom 
(kcal /cm 2 ) 


-0.623 


-0.647 


-0.665 


-0.623 


-0.478 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-1.930 


-1.524 


-1.775 


-2.416 


-2.181 


CM JAN | FEB | MAR | APR , V 




1 AUG | SEP 


, OCT | NOV 








Fig. 12 - Annual equilibrium temperature and thick ness field for Case 4, 
constant salinity (0.09 V.c), F w = 1.5 kcal/cm year, 
bottom temperature - 0.1°C. 



{0.09 ° loo), BOTTOM TEMPERATURE = -0.1 °C. ALL OTHER INPUT DATA AS IN CASE 1. 



-2.30 -0.10 -0.58 -4.87 -10.49 -15.38 -17.91 -11.42 

374.6 355.7 335.5 331.9 330.6 331.4 336.6 348.6 

-0.113 0.013 0.332 0.820 0.797 0.858 0.771 5.965 

-0.285 -0.106 -0.045 -0.029 -0.033 -0.379 -0.577 -4.491 

4.396 4.959 2.710 0.605 0.085 0 0 18.182 

-1.934 -0.903 -1.038 -0.870 -0.797 -1.413 -1.556 -18.338 



Date 



8 June 

28 June 

9 July 
20 August 

29 October 



Snow ablation 
Ice ablation 



Bottom accretion 



Bottom accretion 
Surface ice ablation 



Other Data 

Surface ice ablation' 41.5 cm 

Bottom ablation 4.1 cm 

Bottom accretion 45.6 cm 

Short-wave radiation penetrating 

surface during snow-free period 1.291 kc c - 



Equals net bottom accretion. 
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underside of the ice. Thus their value for F is closer to 6.6 kcal/ 

cm year. As will be seen in Section 4.2.2, an F of this magnitude 

causes the ice to vanish. To test the heat budget of Vowinckel and 

Orvig it was therefore necessary to assume a smaller value for F 

Likewise, it is not possible to estimate from their data how much 

short-wave radiation penetrates to deeper layers of the ice. In the 

following case, F and I were assumed to be as in Case 1 (F = 

2 w ° w 

1.5 kcal/cm year, I q = 17 percent. 

With these values, and the surface heat budget of Vowinckel and 
Orvig, the model predicts 115 cm of surface-ice ablation and complete 
melting of ice during the fifth summer. Comparison of Tables 1 and 
7 shows that Vowinckel and Orvig assume albedos that are up to 20 
percent lower than Marshunova * s ; consequently, the net short-wave 
radiation is 12 kcal/cm 2 year larger. Even increasing I to an un- 
realistically large 50 percent would result in 60 to 70 cm of surface- 
ice ablation and an equilibrium thickness of roughly 1.0 to 1.5 meters. 

The greatest discrepancy between Fletcher's budget and that of 
Vowinckel and Orvig is in the amount of absorbed short-wave radiation. 
As is pointed out in Section 2.4.5, the average heat exchange over a 
large area cannot be used to determine ice growth. If albedos repre- 
sentative of the ice and snow alone are substituted, Vowinckel and 
Orvig' s heat budget should be applicable to the present sea-ice model. 
We therefore repeated the experiment using the albedos suggested by 
Marshunova (Table 1) and the incoming energy fluxes (F^ , F , F , F ) 
given in Table 7. The results of this case (6) are shown in Table 8. 
Case 6 was integrated for a 40-year period, during which the ice grew 
to 540 cm. Although the ice has not strictly achieved equilibrium, 
the initial 1.5 cm difference between surface ablation and net bottom 
accretion would cause the ice to grow only another 20 cm during the 
subsequent 60 years (see Section 4.1.3). The calculated energy fluxes 
and temperatures would not be changed significantly. 

The differences between Case 6 and Case 1 are the result of mark- 
edly different assumptions regarding the turbulent fluxes and incoming 
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Table 7. AVERAGE MONTHLY VALUES FOR THE INPUT DATA AT THE UPPER BOUNDARY 

ACCORDING TO VOWINCKEL AND ORVIG (1966, 1967). 



Input 


Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Year 


Incoming short- 
(kcal/cm 2 ) 

radiation 
(kcal/cm 2 ) 

Implied surface 
albedo 

Incoming long- 
wave radiation 
(kcal/cm 2 ) 

Flux of sensible 
heat 

(kcal/cm 2 ) 
^kcal/cm 2 ) 


0.67 0.73 0.65 0.54 0.43 0.46 0.54 0.63 

0.09 0 0 -0.07 -0.50 -0.44 -0.16 -0.66 -0.69 -0.09 0.09 0 -2.43 



solar energy. Vowinckel and Orvig assume a sensible heat loss at the 

surface 6 kcal/cm^ year greater than Doronin (Table 1) and a latent 
2 

heat loss that is 0.8 kcal/cm year less. Surface ablation in Case 6 
is 75 percent of Case 1 because they assume a smaller and also spec- 
ify smaller turbulent heat losses during the summer. Although decreas- 
ed surface ablation accounts for part of the predicted thickening, the 
greater part is a result of the extremely low winter temperatures at 
the surface. These temperatures are lower than in Case 1 because Vo- 
winckel and Orvig assume that almost no sensible heat reaches the sur- 
face during this time. In July, on the other hand, they assume that 
the heat loss by F g and F^ is less than Doronin 's estimate by 0.8 kcal/ 

cm^ year. During September they estimate that the sum of F and F 

2 s 
is -1.7 kcal/cm , which causes the ice to cool so rapidly that, in 

contrast to observations, the net long-wave radiation balance becomes 

positive . 

The above discussion has shown that none of the values suggested 
by Vowinckel and Orvig for the oceanic heat flux, surface albedo, or 
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Table 8. COMPUTED VALUES FOR CASE 6, INCOMING ENERGY FLUXES AT THE UPPER 
ALL OTHER INPUT DATA AS IN CASE 1. AS DESCRIBED IN 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature ( C) 


-37.95 


-42.49 


-39.45 


-31.24 


-14.45 


Mean ice surface 
temperature ( C) 


-24.53 


-27.44 


-27.78 


-25.06 


-17.60 


Mean ice thickness (cm) 


532.9 


535.5 


539.3 


544.1 


549.3 


Heat flux through surface 
(kcal/cm 2 ) 


0.807 


0.881 


0.659 


0.330 


-0.172 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.245 


-0.335 


-0.407 


-0.449 


-0.453 


Net short-wave radiation 
(kcal /cm 2 ) 


0 


0 


0.368 


1.664 


2.794 


Net long-wave radiation 
(kcal/cm 2 ) 


-1.016 


-0.941 


-1.090 


-1.030 


-1.125 



BOUNDARY FROM VOWINCKEL AND ORVIG (TABLE 7), ALBEDOS FROM MARSHUNOVA (TABLE 1). 
SECTION 4.1.4., THE ICE IS NOT STRICTLY IN EQUILIBRIUM. 



Aug Sep 



-2.75 -0.10 -1.74 -14.21 -20.27 -31.18 -36.06 -22.66 

-9.85 -0.10 -1.41 -10.47 -13.44 -18.10 -22.38 -16.51 

554.0 551.2 538.4 532.8 533.0 532.9 533.2 539.7 

-0.361 -0.103 0.288 0.683 0.550 0.838 0.843 5.244 

-0.403 -0.327 -0.235 -0.162 -0.122 -0.123 -0.164 -3.423 

3.336 4.216 2.555 0.449 0.040 0 0 15.423 

-1.792 -1.081 -0.726 -0.558 -0.499 -0.985 -0.904 -10.633 



Date Onset End 

23 June Snow ablation 

12 July ice ablation 

16 August Surface ice ablation 

10 October Bottom accretion 

20 November Bottom accretion 



Surface ice ablation" 28.5 cm 

Bottom ablation 0.2 cm 

Bottom accretion 30.2 cm 

Short-wave radiation penetrating 2 

surface during snow-free period 0.820 kcal/cm 



Equals net bottom accretion. 
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the turbulent fluxes can be considered representative of conditions 
over the bare ice, although they may be correct in terms of a large- 
scale average. Since ice thickness is not directly a function of en- 
ergy transfer over leads and melt ponds, the heat budget of Vowinckel 
and Orvig cannot be used for realistic calculations of ice growth. 
Likewise, because of a lack of detail in Eq. (34), the sea-ice model 
cannot be used for an explicit evaluation of their assumed energy budget. 
If Vowinckel and Orvig had estimated the separate contributions to 
the 0 term and the effects of leads on the turbulent fluxes, it might 
have been possible to obtain an independent evaluation of their large- 
scale heat budget, assuming that the parameters specified in Case 1 
are characteristic of the bare ice. 

4.2 ANALYSIS OF FURTHER CASES 

The effects of a certain parameter on sea ice may be investigated 
by varying this parameter, while holding all others constant. Since 
the results of Case 1 are compatible with observations, the input data 
for this case have been adopted as the "standard conditions." The 
experiments discussed in this section define how penetrating short- 
wave radiation, oceanic heat flux, and the snow cover influence the 
growth and temperature regime of sea ice. 

A. 2.1 Penetration of short-wave radiation 

The amount of short-wave radiation which penetrates the ice (I q ) 
has been inferred from temperature profiles ( Untersteiner , 1961) , but 
may be in error by as much as a factor of two. Untersteiner (1964) 
has theoretically investigated the effect of 1 q on internal ice tem- 
peratures, but no previous work has related I q to equilibrium thick- 
ness. The experiments described in this section will expand upon 
Untersteiner ' s work and define the role which I q plays in annual mass 
changes . 

It was assumed in Case 1 that 17 percent (about 1.2 kcal/cm 2 year) 
of the net short-wave radiation penetrates the ice during the snow- 
free period. In cases 7 to 10 we have specified that I is 0, 8.5, 
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25.5, and 34 percent of the net short-wave radiation, respectively. 
All other input data are unchanged from Case 1. The results of these 
experiments are given in Figs. 13 to 17 and Tables 9 to 12. 

Contrary to what might be expected, equilibrium thickness actually 
increases as more radiation penetrates the ice (Fig. 17). Because the 
net short-wave radiation remains constant and a greater percentage of 
it is absorbed within the ice, less energy is available for melting 
at the surface. The additional energy absorbed within the ice causes 
slight additional warming in the slab, but it is, for the most part, 
stored as heat of fusion within the diluted brine. Since there is 
little temperature change, the conductive heat flux to the surface 
remains almost unchanged during the ablation season. Tables 9 to 12 
show that, as I q increases from 0 to 34 percent, surface ablation de- 
creases from 52.3 to 24.7 cm/year. Less surface ablation means that the 
ice must thicken to slow bottom accretion and maintain equilibrium. 

Temperatures change surprisingly little as I varies. As I in- 
2 o o 

creases from 0 to 2.5 kcal/cm year, the annual surface temperature 

decreases by only 0.1°C, while the average temperature at the ice/snow 
interface decreases by 1°C. The maximum temperature change within the 
ice is 3°C. Intuitively, it could be argued that an increased energy 
input of 2.5 kcal/cm 2 year into the ice should cause average tempera- 
tures to rise, despite an increase in thickness. An examination of 
Tables 9 to 12 shows that, although downward heat input into the ice 
has increased by 2.5 kcal/cm 2 year, the upward heat flux at the sur- 
face has increased by only 0.4 kcal/cm 2 year. This seeming contradic- 
tion is explained by the 28-cm decrease in bottom accretion, which 
reduces the heat conducted into the ice by about 2.0 kcal/cm 2 year. 
Thus the net energy input is only 20 percent of the increase in I 

o 

and, except during the summer and fall, ice thickness determines the 
temperature profiles. 

The date of the onset of ice and snow ablation does not vary, and 
the duration of surface ablation is nearly the same. As the ice thick- 
ens from 240 to 370 cm, bottom ablation starts about a half -month later 
and is extended by about 2 months. Comparison with cases of equivalent 



- 'Ik - 



Table 9. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature ( C) 


-30.73 


-32.85 


-31.58 


-22.23 


-9.28 


Mean ice surface 
temperature ( C) 


-16.76 


-17.75 


-17.87 


-14.62 


-8.74 




238.3 


246.7 


255.3 


263.4 


269.7 


Heat flux through surface 
(kcal/cm 2 ) 


0.832 


0.876 


0.770 


0.408 




Heat flux through bottom 
(kcal/cm 2 ) 


-0.658 


-0.668 


-0.675 


-0.602 


-0.442 


Net short-wave radiation 
(kcal /cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-2.017 


-1.615 


-1.868 


-2.503 


-2.242 


JAN | FEB | WAR , APR , MA 




, AUG | SEP 


| OCT | NOV 








Fig. 13 — Annual equilibrium temperature and thickness field 
for Case 7, I =0. 



IN CASE 7,1 =0. ALL OTHER INPUT DATA AS IS CASE 1. 



Jun Jul Aug Sep Oct Nov Dec Year 

-0.29 -0.10 -1.14 -10.19 -20.43 -28.71 -30.10 -18.14 

-2.02 -0.10 -0.97 -6.69 -10.91 -14.59 -16.28 ■-10.61 

272.3 250.2 224.5 219.5 219.1 222.9 230.0 242.6 

-0.102 -0.329 . -0.004 0.623 0.768 0.891 0.843 5.595 

-0.263 -0.120 -0.040 -0.013 -0.222 -0.503 -0.641 -4.846 

4.435 4.959 2.789 0.605 0.085 0 0 18.300 

-1.947 -0.903 -0.979 -0.672 -0.767 -1.445 -1.628 -18.587 



Date 



7 June 
28 June 
13 July 
22 August 

7 October 



Snow ablation 
Ice ablation 



Bottom accretion 



Bottom accretion 
Surface ice ablation 



Other Data 

Surface ice ablation 52.3 cm 

Bottom ablation 3.4 cm 

Bottom accretion 55.7 cm 

Short-wave radiation penetrating 

surface during snow-free period 0 

* 

Equals net bottom accretion. 



Table 10. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature ( C) 


-30.85 


-32.98 


-31.71 


-22.33 


-9.35 


Mean ice surface 
temperature ( C) 


-17.31 


-18.32 


-18.44 


-15.15 


-9.16 


Mean ice thickness (cm) 


256 .9 


264.8 


272.9 


280.8 


287.0 


Heat flux through surface 
(kcal/cm 2 ) 


0.806 


0.851 


0.745 


0.384 


0.002 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.629 


-0.638 


-0.648 


-0.592 


-0.447 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-1.992 


-1.590 


-1.843 


-2.480 


-2„224 


1 JAN | FEB | MAR | APR , 




| AUG | SEP 


1 OCT | NOV 








Fig. 



14 — Annual equilibrium temperature and thickness field 
for Case 8, I = 8.5 percent. 



.CASE 8, 1 Q = 8.5%. 



ALL OTHER INPUT DATA AS IN CASE 1. 



Aug Sep 



-2.16 -0.10 -0.97 -6.49 -10.81 -14.76 -16.72 -10.87 

289.9 270.4 247.2 242.6 241.1 243.1 249.2 262.0 

-0.101 -0.158 0.131 0.635 0.772 0.883 0.823 5.773 

-0.276 -0.121 -0.028 0.007 -0.099 -0.408 -0.592 -4.470 

4.433 4.959 2.737 0.605 0.085 0 0 18.246 

-1.945 -0.903 -0.966 -0.685 -0.771 -1.437 -1.608 -18.443 



Date Onset 

7 June Snow ablation 

28 June Ice ablation 
14 July 



Bottom accretio 



21 August Surface ice ablatii 

19 October Bottom accretion 



Surface ice ablation 
Bottom ablation 
Bottom accretion 



Equals net bottom accretion. 



46.5 cm 
4.5 cm 
51.0 cm 



Short-wave radiation penetrating 2 
surface during snow-free period 0.658 kcal/cm 
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Table 11. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


temperature (°C) 


-31.10 


-33.29 


-32.03 


-22.62 


-9.56 


Mean ice surface 
temperature (°C ) 


-18.50 


-19.78 


-20.01 


-16.66 


-10.42 


Mean ice thickness (cm) 


316.4 


322.7 


329.6 


336.5 


342.4 


Heat flux through surface 
(kcal/cm 2 ) 


0.756 


0.789 


0.679 


0.319 


-0.052 


Heat flux through bottom 
(kcal / cm 2 ) 












Net short-wave radiation 
(kcal /cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-1.942 


-1.528 


-1.777 


-2.414 


-2.169 


CM 


A¥ ,,l JUS 1 m 


„l,., Aur ' 1 S!P 


1 OCT j NOV 


DEC 
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Fig. 15 -- Annual equilibrium temperature and thickness field 
for Case 9, I = 25.5 percent. 
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IN CASE 9,1 = 25.5%. ALL OTHER INPUT DATA AS IN CASE 1. 



Sep 



-0.35 -0.10 -0.95 -9.80 -20.31 -28.72 -30.32 -18.26 

-^.63 -0.10 -0.61 -5.71 -10.50 -14.75 -17.35 -11.42 

345.8 332.4 316.5 313.6 311.6 310.1 311.6 324.1 

-0.130 0.014 0.325 0.726 0.794 0.888 0.798 5.906 

-0.308 -0.159 -0.056 -0.012 0.005 -0.090 -0.351 -3.596 

4.359 4.959 2.710 0.605 0.085 0 0 18.145 

-1.929 -0.903 -1.030 -0.776 -0.794 -1.442 -1.583 -18.288 



Date 



8 June 
29 June 
22 July 
17 August 
21 November 



Snow ablation 
Ice ablation 



Bottom accretion 



Bottom accretion 
Surface ice ablation 



Surface ice ablation 32.8 cm 

Bottom ablation 5.8 cm 

Bottom accretion 38.6 cm 
Short-wave radiation penetrating 

surface during snow-free period 1.910 kcal/cm 2 



Equals net bottom accreti< 



Table 12. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature ( C) 


-31.15 


-33.41 


-32.20 


-22.79 


-9.68 


Mean ice surface 
temperature ( C) 


-18.72 


-20.33 


-20.76 


-17.48 


-11.19 


Mean ice thickness (cm) 


359.8 


364.3 


370.0 


376.1 


381.7 


Heat flux through surface 
(kcal/cm 2 ) 


0.746 


0.765 


0.647 


0.282 


-0.085 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.360 


-0.462 


-0.513 


-0.516 


-0.438 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-1.932 


-1.504 


-1.744 


-2.378 


-2.136 


JAN 1 ,LB , i MAR 1 APR !. MA 


Y | JUN | JUL 


,1 AUG 1 SEP .. 


OCT | NOV 
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Fig. 16 -- Annual equilibrium temperature and thickness field 
for Case 10, I =34 percent. 
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IN CASE 10, I = 34%. ALL OTHER INPUT DATA AS IN CASE 1. 



Aug Sep 



-2.92 -0.10 -0.50 -4.45 -9.95 -14.50 -17.26 -11.51 

385.4 375.6 364.1 362.3 360.6 358.8 357.8 368.1 

-0.141 0.037 0.355 0.870 0.829 0.901 0.802 6.006 

-0.323 -0.187 -0.082 -0.035 -0.010 -0.009 -0.144 -3.080 

4.354 4.959 2.710 0.605 0.085 0 0 18.140 

-1.922 -0.903 -1.043 -0.920 -0.828 -1.455 -1.588 -18.354 



Date 



8 June 
29 June 

1 August 
14 August 
14 December 



Snow ablation 
Ice ablation 



Bottom accretion 



Bottom accretion 
Surface ice ablation 



Other Data 

Surface ice ablation 24.7 cm 

Bottom ablation 6.0 cm 

Bottom accretion 30.7 cm 

Short-wave radiation penetrating 

surface during snow-free period 2.547 kcal/cm^ 



Equals net bottom accretion. 
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thicknesses (e.g., Case 13) shows that bottom accretion is initiated 
approximately a month later for a large I q . 

As was pointed out in Section 4.1.1, the model predicts more rapid 
cooling in the fall than observations indicate. One possible reason 
for this is that we have underestimated I . Comparison of Figs. 6, 
7, and 16 shows that, when I q is 34 percent (Case 10), the fall tem- 
peratures are in much better agreement with empirical data than in 
Case 1; however, ice ablation at the surface is too small. 

It is possible to increase surface ablation, while maintaining 
a large I q simply by decreasing the albedo of the melting ice (o^) ; 
a ± is sufficiently in doubt that a small percentage change is accept- 
able. In Case 11, a has been reduced from 0.64 to 0.58 and I spec- 

1 o r 

ified as 34 percent. Results are given in Table 13 and Fig. 18. A 

comparison of Cases 1 and 11 shows nearly the same surface ablation, 
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Fig. 17 



— Average equilibrium thickness as a function of the per- 
centage of net short-wave radiation which penetrates the 
ice during the snow-free period. 



-85- 



but solar energy absorbed within the ice is greater by 1.9 kcal/cm yr. 
Since mass changes at the surface have been separated from I q , it now 
becomes possible to define how the ice responds to I q alone. Even 
though surface ablation remains nearly constant, bottom ablation in 
Case 11 is more than double that of Case 1. To maintain equilibrium, 
bottom accretion must therefore increase, but this may occur only if 
the ice thins. Average thickness in Case 11 is 229 cm, that is, 60 cm 
less thickness for an annual energy greater by only 1.9 kcal/cm . Fig- 
ure 23 shows that an equivalent increase in F w would cause the equilib- 
rium thickness to decrease from 288 cm to 140 cm, illustrating the 
greater efficiency of heat addition from the ocean. 

In Case 11 the rate of cooling within the ice in autumn is con- 
siderably less than in Case 1. Surface temperatures are higher in the 
winter, but this is due to thinner ice. Fall temperature patterns in 
Fig. 18 are similar to the observed ones (Fig. 7); however, the cooling 
rate in the upper meter of the ice is still too high September through 
December. It does not appear that a further increase in I would be 
able to duplicate the pattern of Fig. 7. Because the greatest dis- 
crepancy now occurs in the upper portion of the ice, the most likely 
explanation for the slower cooling in autumn is the latent-heat re- 
lease from freezing melt ponds, possibly combined with an underestima- 
tion of I . 

Although Untersteiner (1961) calculated an I q of 34 percent dur- 
ing the height of the melt season and Chernigovskii (1963) appears to 
have measured an I of 47 percent, there are several objections to 
values this large. When I q is large, the combination of a thickness 
of 300 cm and a surface ablation of 40 cm is impossible. If we accept 
less ice ablation, it is possible to have a reasonable ice thickness 
and a large I q . Untersteiner, for example, measured only 32 cm of 
ice melt, which is consistent with his ice thickness and I measure- 
ments. Ablation during the following year was 53 cm, but no I q calcu- 
lations were reported. 

These considerations demonstrate the need for careful studies of 
the relationship between I , net short-wave radiation, and surface ab- 
lation. Observations of melt ponds after the fall freeze-up would 
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Table 13. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 11, 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature ( C) 


-30.58 


-32.71 


-31.20 


-22.13 


-9.22 


temperature (°C) 


-16.18 


-17.23 


-17.37 


-14.15 


-8.36 


Mean ice thickness (cm) 


221.8 


230.6 


239.6 


248.2 


254.6 


Heat flux through surface 
(kcal/cm 2 ) 


0.863 


0.905 


0.796 


0.431 


0.036 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.683 


-0.696 


-0.703 


-0.622 


-0.440 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-2.048 


-1.643 


-1.894 


-2.526 


-2.258 




Fig. 18 -- Annual equilibrium temperature and thickness field f< 
Case 11, I =34 percent, surface ice albedo = 0.58. 
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I = 34%, ICE ALBEDO = 0.58. ALL OTHER INPUT DATA AS IN CASE 1. 



Jul 



-0.30 -0.10 -0.86 -9.21 -19.80 -28.30 -29.90 -17.88 

-'.84 -0.10 -0.48 -4.26 -8.53 -12.92 -15.47 -9.74 

257.4 236.8 217.2 212.6 209.1 207.6 213.2 229.1 

-0.091 0.040 0.364 0.881 0.914 0.977 0.886 7.003 

-0.250 -0.034 0.082 0.100 0.096 -0.249 -0.632 -4.032 

4.817 5.785 3.098 0.605 0.085 0 0 19.818 

-1.944 -0.903 -1.057 ' -0.931 -0.914 -1.532 -1.671 -19.322 



Date 



7 June 
26 June 

3 July 
18 August 
10 November 



Snow ablation 
Ice ablation 



Bottom accretion 



Bottom accretion 
Surface ice ablation 



Other Data 

Surface ice ablation 39.6 cm 

Bottom ablation 12.1 cm 

Bottom accretion 51.7 cm 

Short-wave radiation penetrating 

surface during snow-free period 3.176 kcal/cm 2 



Equals net bottom accretion. 
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also be useful. Because it is still impossible to say which combina- 
tion of I q and surface ablation is truly representative, the succeed- 
ing cases were integrated with an I q of 17 percent, conceding, however, 
that I may be larger. 

4.2.2 Oceanic heat flux 

Since the pioneering work of Malmgren (1927), the turbulent trans- 
fer of heat from the ocean to the ice has been subject to numerous es- 
timates and speculations (see Chapters I and II) . Wide variations in 
F w are an integral part of several theories of the ice ages, and num- 
erous schemes for artificially influencing the extent of the pack ice 
depend upon changing F w ( Toporkov , 1963; Fletcher , 1965). The feasi- 
bility of any of these proposals is difficult to establish because of 
a lack of quantitative knowledge regarding the sensitivity of ice re- 
sponse to F w> It is not known how large the heat flux must become 
before the ice would vanish, and disagreement still exists as to the 
present magnitude of F^. In this section an attempt will be made first 
to describe the total response of the ice to changes in F w , then to 
determine an upper limit for F w , and finally to define a range of con- 
ditions under which the ice could currently exist. 

Five cases were integrated using annual values for F w of 0, 0.75, 
3.0, 4.5, and 6.0 kcal/cm 2 year (Cases 12 to 16, respectively). Be- 
cause no obvious inconsistencies had resulted from assuming that F 
remained constant throughout the year in Case 1, the other cases were 
computed under the same assumption. All other fluxes were left un- 
changed. The temperature fields and monthly averages for Cases 12 to 

15 appear in Figs. 19 to 22 and Tables 14 to 17. When F =6.0 kcal/ 

2 w 
cm year, the ice vanishes before equilibrium is established and con- 
sequently no other results are presented for this case. 

Figure 23 shows a graph of equilibrium thickness vs . F ; curves 

w 

for seasonal maximum and minimum thickness and average thickness are 
drawn. When the ice is thicker, the average value is closer to the 
minimum, reflecting the decrease in bottom ablation. Under present 
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surface conditions, the theoretical thermodynamic limit of ice thick- 
ness is about 560 cm, whereas the ice should vanish if F increases 

2 w 
to five or six kcal/cm year. 

An examination of surface conditions shows that variations in 
oceanic heat flux have relatively little direct influence on processes 
occurring at the upper boundary. The surface temperature of ice 90 cm 
thick is up to 12°C higher than the surface temperature of 560-cm ice, 
but because the snow is an effective thermal insulator, average snow 
surface temperatures are only two or three degrees higher. Likewise, 
because of the snow cover, the onset of snow and ice ablation changes 
by only three days. A somewhat surprising result is that the upper 
surface ice ablation varies by only a few centimeters. Tables 2 and 
14 to 17 show that the quantity of heat conducted to the surface dur- 
ing the snow-free period is similar in all cases, indicating that tem- 
perature gradients during the summer are generally determined by the 
penetrating short-wave radiation, rather than by the ice thickness. 

Figure 24 shows cumulative bottom accretion vs. time and illus- 
trates the earlier onset of, and greater amount of, bottom ablation 
associated with higher F^. Bottom accretion, unexpectedly, also be- 
gins earlier in the fall when F w is larger. Examination of the tem- 
perature fields verifies that, in thin ice, the lower portion responds 
more rapidly to temperature changes at the surface. Thus the steeper 
temperature gradients within the ice outweigh the effects of a higher 
F w and result in earlier bottom accretion. 

Most field observations indicate an average ice thickness between 
250 and 350 cm, which corresponds in Fig. 23 to values of F between 
1.0 and 2.0 kcal/cm year. These values refer to an "effective" F . 
It is probable that more heat is carried into the basin but is lost 
through leads. Because uncertainties that could also affect thickness 
exist in the other input data, it is not possible to assign a specific 
value to F w - However, on the basis of the temperature results, it is 
unlikely that errors in the energy fluxes at the upper surface could 
modify the thickness results by more than ± 0.5 m. This being the 
case, the estimates of oceanic heat flux by Crary, Badgley, and Unter- 
steiner must be very close to reality and the earlier ones considerably 
in error. 



Table 14. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature (°C) 


-31.64 


-33.90 


-32.70 


-23.27 


-10.09 


Mean ice surface 
temperature ( C) 


-20.88 


-22.47 


-23.05 


-19.94 


-13.67 


Mean ice thickness (cm) 


557.0 


560.0 


564.0 


568.8 


574.0 


Heat flux through surface 
(kcal/cm 2 ) 


0.646 


0.669 


0.545 


0.173 


-0.191 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.156 


-0.225 


-0.283 


-0.326 


-0.334 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-1.832 


-1.408 


-1.643 


-2.268 


-2.030 




Fig. 19 -- Annual equilibrium temperature and thickness field 
for Case 12, F =0. 



CASE 12, F w = 0. 



ALL OTHER INPUT DATA AS IN CASE 1. 



Aug 



-4.33 
579.0 



4.237 
-1.895 



-0.10 
565.8 



-0.196 -0.072 
-0.309 -0.252 



4.959 
-0.904 



-1.09 -10.47 -20.86 



-0.81 
551.1 

0.266 
-0.185 

2.710 
-0.994 



-7.31 
550.6 

0.553 
-0.131 

0.605 
-0.602 



-12.58 
552.4 

0.668 
-0.097 

0.085 
-0.668 



-29.26 -30.85 



-17.08 
553.8 



-19.64 
555.1 



0.776 0.689 
-0.083 -0.100 



-13.49 
561.0 

4.525 
-2.478 
18.023 
-17.048 



9 June Snow ablation 

1 July Ice ablation 



Bottom accretion 



Bottom accretion 
Surface ice ablation 



Other Data 

* 

Surface ice ablation 38.7 cm 

Bottom ablation 0 

Bottom accretion 38.7 cm 

Short-wave radiation penetrating 

surface during snow-free period 1.222 kcal/cm 2 



Equals net bottom accretion. 



Table 15. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 13, 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


temperature (°C) 


-31.33 


-33.55 


-32.32 


-22.89 


-9.77 


temperature (°C ) 


-19.50 


-20.95 


-21.31 


-18.01 


-11.68 


Mean ice thickness (cm) 


384.5 


390.0 


396.4 


403.1 


409.3 


Heat flux through surface 
(kcal/cm 2 ) 


0.709 


0.737 


0.622 


0.259 


-0.106 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.380 


-0.449 


-0.489 


-0.486 


-0.430 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-1.895 


-1.476 


-1.720 


-2.354 


-2 0 115 


JAN | FIB | MAR | APR , MA 




1 AUG 1 SEP 


, OCT , NOV 








Fig. 20 — Annual equilibrium temperature and thickness field for 
Case 13, F =0.75 kcal/cm year. 
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F =0.75 kcal/cm" year. ALL OTHER INPUT DATA AS IN CASE 1. 



Jun 


Jul 


Aug 


Sep 


Oct 


Nov 


Dec 


Year 


-0.39 


-0.10 


-1.08 


-10.23 


-20.59 


-28.97 


-30.55 


-18.48 


-3.20 


.-0.10 


-0.80 


-6.72 


-11.56 


-15.81 


-18.31 


-12.33 


413 .8 
















-0.156 


-£.042 


0.270 


0.614 


0.730 


0.837 


0.752 


5.225 


-0.325 


-0.208 


-0.113 


-0.061 


-0.039 


-0.084 


-0.237 


-3.302 


4.319 


4.959 


2.710 


0.605 


0.085 


0 


0 


18.105 


-1.918 


-0.903 


-0.996 


-0.663 


-0.730 


-1.391 


-1.538 


-17.699 



Date 



8 June 
29 June 
19 August 
13 September 
10 November 



Snow ablation 
Ice ablation 



Bottom accretion 



Surface ice ablati* 
Bottom accretion 



Other Data 

Surface ice ablation 39.9 cm 

Bottom ablation 0.5 cm 

Bottom accretion 40.4 cm 
Short-wave radiation penetrating 

surface during snow-free period 1.256 kc; 



Equals net bottom accretion. 



Table 16. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 14, 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature (°C) 


-30.03 


-32.15 


-30.90 


-21.68 


-8.95 


Mean ice surface 
temperature (°C) 


-13.73 


-14.73 


-14.86 


-11.90 


-6.74 


Mean ice thickness (cm) 


159.3 


168.0 


177.0 


184.8 


189.5 


Heat flux through surface 
(kcal/cm 2 ) 


0.976 


1.017 


0.908 


0.532 


0.107 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.808 


-0.821 


-0.814 


-0.670 


-0.425 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation ] 
(kcal/cm 2 ) j 


-2.162 


-1.756 


-2.006 


-2.627 


-2.329 



CM 




Fig. 21 -- Annual equilibrium temperature and thickness field for 
Case 14, F =3.0 kcal/cm 2 year. 



F =3.0 kcal/cm^ year. ALL OTHER INPUT DATA AS IN CASE 1. 



-0.26 -0.10 -1.11 -9.69 -20.06 -28.15 -29.41 -17.71 

-1.37 -0.10 -0.83 -5.39 -9.50 -12.23 -13.32 -8.72 

189.4 168.1 144.2 136.2 135.3 141.6 150.3 161.9 

-0.079 -0.098 0.206 0.753 0.853 1.010 0.988 7.172 

-0.194 0.027 0.145 0.093 -0.504 -0.766 -0.830 -5.566 

4.469 4.959 2.710 0.605 0.085 0 0 18.255 

-1.956 -0.903 -0.986 -0.803 -0.852 -1.563 -1.772 -19.714 



Date 

7 June 

7 June 
28 June 
18 August 

3 October 



Onset 
Snow ablation 

Ice ablation 

Bottom accretion 



Bottom accretion 



Surface ice ablatii 



Surface ice ablation 40.1 cm 

Bottom ablation 16.7 cm 

Bottom accretion 56.7 cm 

Short-wave radiation penetrating 

surface during snow-free period 1.308 kcal/cm 2 



Equals net bottom accretioi 



Table 17. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 15, 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature ( C) 


-29.21 


-31.35 


-30.15 


-21.12 


-8.66 


Mean ice surface 
temperature (°C ) 


-10.11 


-11.17 


-11.42 


-9.07 


-5.01 


Mean ice thickness (cm) 


93.8 


103.3 


112.8 


120.2 


123.0 


Heat flux through surface 
(kcal/cm 2 ) 


1.144 


1.178 


1.061 


0.658 


0.182 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.976 


-0.993 


-0.952 


-0.724 


-0.389 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-2.330 


-1.917 


-2.159 


-2.755 


-2.405 




Fig. 22 -- Annual equilibrium temperature and thickness field for 
Case 15, F w = 4.5 kcal/cm year. 



=4.5 keal/cm year. ALL OTHER INPUT DATA AS IN CASE 1. 



-9.58 -19.39 



-27.31 -28.54 



-0.98 
120.0 

-0.069 
-0.095 
4.503 
-1.964 



-0.10 
95.0 

-0.195 
0.218 
4.959 

-0.903 



-0.80 
67.2 

0.134 
0.381 
2.710 
-0.991 



-5.05 
60.0 

0.781 
-0.517 

0.605 
-0.830 



-6.89 
65.4 

1.008 
-0.858 

0.085 
-1.007 



-8.54 
74.0 

1.189 
-0.999 



-9.49 
84.1 

1.169 
-1.015 



-6.55 
93.2 

8.241 
-6.919 
18.289 
-20.957 



Date 

16 May 

6 June 
27 June 

17 August 

7 September 



Snow ablation 
Ice ablation 



Bottom accretion 



End 

Bottom accretion 



Surface ice ablation 



Surface ice ablation 37.8 cm 

Bottom ablation 25.5 cm 

Bottom accretion 63.3 cm 

Short-wave radiation penetrating 

surface during snow-free period 1.326 kc< 



Equals net bottom accretion. 




Fig. 24 -- Net bottom accretion as a function of time for equilibrium 

ice with various oceanic heat fluxes (F ) . 

vr 
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4.2.3 Snow cover 

A snow cover insulates the ice from the cold air and reflects 
much of the incoming short-wave radiation. The effects of snow on 
surface ablation and ice temperatures suggest that small changes in 
snow thickness might result in significant changes in equilibrium 
thickness. The series of experiments described in this section indi- 
cates that this is not the case. 

It was assumed in Case 1 that the maximum snow depth under pre- 
sent conditions is 40 cm. To investigate the effects of variable 
snow thickness, six cases were integrated specifying maximum snow 
depths of 0, 20, 60, 80, 100, and 120 cm. Snow was added in incre- 
ments as described in 2.4.4, except that each accumulation increment 
was multiplied by a constant factor (0.5, 1.5, etc.) to bring the 
final snow depth up to the desired total. Except for the no-snow 
Case (17), all assumptions were consistent with Case 1. When the 
snow was absent during the entire year, it was assumed that = 0.75 
if the surface temperature was below freezing and 0.64 during periods 
of surface melting. Figures 25 to 30 and Tables 18 to 23 present the 
results of these experiments . 

Figure 30 shows the average equilibrium thickness as a function 
of maximum annual snow depth. It is surprising that, for snow depths 
between 0 and 70 cm, the average ice thickness remains almost un- 
changed; above 70 cm, the thickness increases rapidly. To explain 
the unusual shape of this curve, it is necessary to examine the three 
major variables which are influenced by the amount of snow: ice tem- 
perature, surface ice ablation, and I q . A thinner snow cover allows 
greater cooling of the ice during the winter, promoting bottom accre- 
tion. If ice ablation at the surface remained constant, greater bot- 
tom accretion would lead to thicker equilibrium ice. However, when 
the maximum snow depth is less, the snow cover is removed earlier, 
thereby decreasing the average albedo and prolonging the period of 
ice ablation. Also, less energy is required to remove the thinner 
snow cover. With these effects ice ablation is greater and more 
short-wave radiation is absorbed within the ice. Greater ice ablation 
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at the surface tends to thin the ice but, as we saw in Section 4.2.1, 
I o interacts witn the ice in a more complex manner. Greater I causes 
more bottom ablation, tending to thin the ice, while simultaneously 
depriving the surface of energy which could otherwise be used for melt- 
ing. Thus there are four opposing effects: two which promote thick- 
ening and two which promote thinning. Figure 30 results from the bal- 
ance of these effects. When maximum snow depth lies between 0 and 
70 cm, there is a mutual compensation between the thickening effects 
and the thinning effects; when the annual snow depth exceeds 70 cm, 
decreasing ice ablation and a smaller I dominate the balance and the 
ice thickens to maintain the balance. 

These points are illustrated by the contrast between Case 1 and 
Case 17, the no-snow case. Ice surface temperatures in Case 17 are 
consistently lower, the maximum difference being about 10°C, while 
the annual average differs by 5°C. As may be seen from Fig. 25, tem- 
perature gradients within the ice are very steep during the winter 
months, resulting in 30 cm more bottom accretion than in Case 1. Be- 
cause there is no snow to melt and because a. is smaller than a , sur- 

i s 
face ice ablation is greater by 28 cm. Although net short-wave radia- 
tion is 4.5 kcal/cm 2 year higher, 2.6 kcal/cm 2 year of this excess is 
absorbed within the ice, producing slightly more bottom ablation. Ice 
thickness during the spring is somewhat greater than in Case 1, but in 
the fall and early winter it is appreciably less. Thus, despite large 
differences in the mass changes at the boundaries, the annual energy 
balance and thickness of the slab are only slightly affected. A sim- 
ilar analysis of Cases 20 to 22 shows how less net short-wave radia- 
tion and the greater snow cover combine to check surface-ice ablation 
and I q . The smaller energy input to the slab rapidly becomes the dom- 
inant influence and results in the predicted thickening of the ice. 

Under the assumed energy input, the maximum snow depth which may 
occur, without producing an annual surplus, is 120 cm. If snow accu- 
mulation exceeds this amount, equilibrium ice will grow from above 
and ablate from below. The metamorphism of snow into ice cannot be 
handled by the present model. From this series of experiments, it 
appears that present variations in annual snowfall over the Central 
Arctic have little effect on average ice thickness. 
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4.3 EVALUATION OF METHODS TO INFLUENCE 
THE ICE COVER ARTIFICIALLY 

Man may be able to induce large-scale changes in climate by arti- 
ficially altering the horizontal extent of pack ice in the Arctic. Al- 
though the exact climatic consequences of such changes are unknown, 
quantitative results from the sea-ice model may be used to evaluate 
the feasibility of various schemes to accomplish these changes. 

Barring changes in the energy balance at the surface, the effec- 
tive oceanic heat flux must be increased by nearly 400 percent (see 
Section 4.2.2) if the ice in the Central Arctic is to vanish. It has 
been proposed that, by pumping water out of the Arctic Ocean across a 
dam at Bering Straits, the increased inflow of Atlantic water could 
dissipate the ice cover ( Toporkov , 1963). As Badgley (1961) has point- 
ed out, even doubling the heat flux by this method is beyond current 
technological capabilities. Another scheme suggests covering the Nor- 
wegian Sea with chemicals to retard cooling, thereby producing a tem- 
perature increase of 3 to 4°C in the entering Atlantic waters ( Maka- 
rova et al. , 1962). This might remove ice in the peripheral seas, at 
least, but as the area of open water increases, more and more heat 
will be lost. Unless chemicals are spread across these areas also, 
it seems unlikely that this method can seriously affect the central 
basin. The conclusion must be that it is not currently feasible to 
remove the pack ice by artificially changing F . It is beyond the 
scope of this paper to try to estimate the necessary rise in sea level 
required for a four-fold increase in F w ; this estimate should be made, 
however, since the ice-age theories of Ewing and Donn (1956, 1958, 
1966) rest on the assumption that the pack ice vanished in response 
to a natural rise in sea level. 

It is clear from Cases 17 to 22 that artificial suppression of 
snowfall would be an ineffective means of decreasing the ice pack be- 
cause of the compensating increase in winter bottom accretion. Early 
removal of the snow, perhaps by annual flooding of the floes with sea 
water in May or early June, offers the only possibility of using the 
snow to decrease ice thickness . 



Table 18. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN 
SURFACE TEMPERATURE IS LESS THAN 272.9 °K. 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature ( C) 












Mean ice surface 
temperature ( C) 


-28.86 


-30.80 


-29.79 


-21.01 


-7.48 


Mean ice thickness (cm) 


300.2 


312.7 


325.4 


337.1 


345.2 


Heat flux through surface 
(kcal/cm 2 ) 


1.220 


1.291 


1.076 


0.667 


0.301 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.930 


-0.936 


-0.930 


-0.783 


-0.496 


Net short-wave radiation 
(kcal /cm 2 ) 


0 


0 


0.597 


2.490 


4.346 


Net long-wave radiation 
(kcal/cm 2 ) 


-2.404 


-2.029 


-2.235 


-2.781 


-2.722 




Fig. 25 



- Annual equilibrium temperature and thickness field for 
Case 17, maximum snow depth =0. 



CASE 17, MAXIMUM SNOW DEPTH =0, ICE ALBEDO = 0.75 IF THE 
ALL OTHER INPUT DATA AS IN CASE 1. 



-0.22 
337.3 

-0.027 
-0.232 
6.313 
-1.954 



-0.10 
301.9 

0.008 
-0.084 

4.959 
-0.903 



-0.72 
280.9 

0.305 
-0.013 

2.909 
-1.098 



-7.89 
277.3 



0.945 
-1.279 



-17.52 
275.0 



1.134 1.438 
0.020 -0.013 



-25.86 
277.2 

1.495 
-0.583 



0.142 0 
-1.451 -2.054 



-28.15 
287.7 

1.246 
-0.910 



-16.53 
304.8 

10.154 
-5.889 

22.700 
-22.946 



Snow ablation 
6 June Ice ablation 

3 July Bottom accretion 

19 August Surface ice ablation 

28 October Bottom accretion 



Other Data 

* 

Surface ice ablation 68.5 cm 

Bottom ablation 6.4 cm 

Bottom accretion 74.9 cm 

Short-wave radiation penetrating ^ 
surface during snow-free period 3.859 kcal/cm 2 



Equals net bottom accretion. 



Table 19. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 18, 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature ( C) 


-30.32 


-32.39 


-31.34 


-22.51 


-9.84 


Mean ice surface 
temperature ( C) 


-22.67 


-24.08 


-24.10 


-19.33 


-11.04 


Mean ice thickness (cm) 


312.3 


320.9 


329.8 


338.6 


345.9 


Heat flux through surface 
(kcal/cm 2 ) 


0.916 


0.969 


0.821 


0.347 


-0.125 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.663 


-0.684 


-0.701 


-0.657 


-0.515 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.892 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-2.101 


-1.708 


-1.917 


-2.439 


-2.095 



MAXIMUM SNOW DEPTH = 20 cm. ALL OTHER INPUT DATA AS IN CASE 1. 



Jun Jul Aug Sep Oct Nov Dec Year 









-9 .65 




-27 . 78 


-29.57 




















348.9 


324.0 


303.8 


300.8 


299.0 


299.2 


304.5 


319.0 


-0.225 


-0.007 


0.287 


0.763 


0.966 


1.087 


0.954 


6.753 


-0.334 


-0.164 


-0.060 


-0.016 


-0.028 


-0.298 


-0.585 


-4.705 


4.851 


4.959 


2.710 


0.605 


0.085 


0 


0 


18.637 


-1.886 


-0.903 


-1.006 


-0.813 


-0.966 


-1.642 


-1.738 


-19.214 



Date Onset End 

9 June Snow ablation 

21 June Ice ablation 

23 July Bottom accretion 

20 August Surface ice ablation 

2 November Bottom accretion 



Other Data 

* 

Surface ice ablation 50.1 cm 

Bottom ablation 4.3 C m 

Bottom accretion 54.4 C m 

Short-wave radiation penetrating 

surface during snow-free period 1.557 kcal/cm 2 



Equals net bottom accretion. 



Table 20. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 19, 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 












temperature (°C) 


-31.57 


-33.77 


-32.41 


-22.70 


-9.39 


Mean ice surface 












temperature ( C) 


-15.06 


-15.98 


-16.19 


-13.69 


-8 88 


Mean ice thickness (cm) 


276.8 


282.6 


288.7 


294.9 


300.0 


Heat flux through surface 












(kcal/cm 2 ) 


0.661 


0.694 


0.603 


0.302 


-0.007 


Heat flux through bottom 












(kcal/cm 2 ) 


-0.483 


-0.505 


-0.528 


-0.498 


-0.399 


Net short-wave radiation 












(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 












(kcal/cm 2 ) 


-1.846 


-1.433 


-1.701 


-2.398 


-2.215 




Fig. 26 -- Annual equilibrium temperature and thickness field for 
Case 19, maximum snow depth = 60 cm. 



MAXIMUM SNOW DEPTH = 60 cm. ALL OTHER INPUT DATA AS IN CASE 1. 



Jun Jul Aug Sep Oct Nov Dec Year 



-0.29 


-0.10 


-1.18 


-10.39 


-20.94 


-29.46 


-30.89 


-18.59 


-2.42 


-0.10 


-0.81 


-5.51 


-8.97 


-12.17 


-14.24 


-9.50 


303.2 


294.1 


277.9 


271.7 


269.6 


269.0' 


271.7 


283.3 


-0.081 


-0.052 


0.226 


0.572 


0.650 


0.735 


0.682 


4.984 


-0.267 


-0.133 


-0.037 


0.004 


-0.008 


-0.185 


-0.391 


-3.429 


4.087 


4.839 


2.710 


0.605 


0.085 


0 


0 


17.752 


-1.948 


-0.909 


-0.967 


-0.621 


-0.649 


-1.290 


-1.467 


-17.455 



Date Onset End 

7 June Snow ablation 

6 July Ice ablation 
17 July Bottom accretion 

19 August Surface ice ablation 

9 November Bottom accretion 



Other Data 

Surface ice ablation" 30.2 cm 

Bottom ablation 5.5 cm 

Bottom accretion 35 # 7 cm 

Short-wave radiation penetrating 

surface during snow- free period 1.042 kcal/cm 2 



Equals net bottom accretion. 



Table 21. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 20, 



Variable 






Mar 


Apr 


May 


Mean snow surface 












temperature ( C) 


-32.08 


-34.34 


-32.93 


-22.96 


-9.47 


Mean ice surface 












temperature ( C) 


-13.55 


-14.51 


-14.84 


-12.91 


-8.94 


Mean ice thickness (cm) 


310.3 


314.1 


318.6 


323.3 


327.4 


Heat flux through surface 












(kcal/cm 2 ) 








0„247 


-0.028 


Heat flux through bottom 












(kcal/cm 2 ) 


-0.344 


-0.393 


-0.425 


-0.416 


-0.360 


Net short-wave radiation 












(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 












(kcal/cm 2 ) 


-1.744 


-1.322 


-1.599 


-2.341 


-2.194 


JAN | FEB | MAR | APR | M 


*Y | JUN , JUL 


| AUG | SEP 


| OCT | NOV 


| DEC 






Fig. 



27 -- Annual equilibrium temperature and thickness field 
for Case 20, maximum snow depth = 80 cm. 



MAXIMUM SNOW DEPTH = 80 cm. ALL OTHER INPUT DATA AS IN CASE 1. 

Jim Jul Aug Sep Oct Nov Dec Year 

-0.28 -0.10 -1.23 -10.71 -21.40 -30.00 -31.41 -18.91 

-2.92 -0.10 -0.77 -5.19 -8.12 -10.72 -12.62 -8.77 

330.4 328.0 316.7 310.0 308.3 306.9 307.5 316.8 

-0.083 -0.040 0.213 0.489 0.543 0.624 0.577 4.186 

-0.265 -0.157 -0.069 -0.025 -0.013 -0.085 -0.238 -2.790 

3.940 4.516 2.710 0.605 0.085 0 0 17.282 

-1.949 -0.918 -0.954 -0.538 -0.543 -1.177 -1.362 -16.642 



Date 



7 June 
15 July 
24 July 
19 August 
24 November 



Snow ablation 
Ice ablation 



Bottom accretion 



Bottom accretion 
Surface ice ablatio: 



Other Data 

Surface ice ablation" 20.2 < 
Bottom ablation 4.8 c 

Bottom accretion 25.0 < 
Short-wave radiation penetrating 

surface during snow-free period 0.775 



Equals net bottom accretion. 



Table 22. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 21, 



Variable 










May 


Mean snow surface 
temperature ( C) 


-32.51 


-34.81 


-33.40 


-23.29 


-9.60 


Mean ice surface 
■ temperature ( C ) 






-14.61 


-13.27 


-9.96 


Mean ice thickness (cm) 


406.2 


407.5 


409.7 


412.4 


415.3 


Heat flux through surface 
(kcal/cm 2 ) 


0.473 


0.492 


0.409 


0.174 


-0.061 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.176 


-0.241 


-0.285 


-0.308 


-0.303 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-1.658 


-1.231 


-1.506 


-2.267 


-2.160 


CM 

JAN | FEB j MAR | APR , MAY 


1 JUN , JUL 


1 AUG | SEP 


OCT | NOV 


DEC 






Fig. 28 -- Annual equilibrium temperature and thickness field 
for Case 21, maximum snow depth = 100 cm. 



JIAXIMUM SNOW DEPTH = 100 cm. ALL OTHER INPUT DATA AS IN CASE 1. 



Sep 



-3.82 -0.10 -0.76 -5.27 -8.06 -10.36 -12.07 -8.78 

417.8 419.0 415.8 409.0 408.0 406.8 406.1 411.1 

-0.094 -0.023 0.188 0.410 0.453 0.524 0.487 3.432 

-0.262 -0.194 -0.124 -0.077 -0.052 -0.056 -0.101 -2.179 

3.874 4.140 2.710 0.605 0.085 0 0 16.841 

-1.946 -0.928 -0.940 -0.459 -0.453 -1.079 -1.271 -15.898 



Date 



7 June 

25 July 
14 August 
19 August 

26 December 



Snow ablation 
Ice ablation 



Bottom accretion 



Bottom accretion 
Surface ice ablation 



Surface ice ablation 10.6 cm 
Bottom ablation 3^4 cm 

Bottom accretion ^ q cm 
Short-wave radiation penetrating 

surface during snow-free period 0.502 kcal/ci 



Equals net bottom accretion. 



Table 23. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 22, 



Variable 


Jan 


Feb 


Mar 


Apr 




Mean snow surface 
temperature (°C) 


-32.97 


-35 .29 






■ 

-9.70 


Mean ice surface 
temperature ( C) 


-14.08 








-11.91 


Mean ice thickness (cm) 


702.0 


701.5 


701.2 


701.2 


701.4 


Heat flux through surface 
(kcal/cm2) 


0.381 


0.398 


0.323 


0.108 


-0.085 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.092 


-0.099 


-0.113 


-0.131 


-0.151 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-1.567 


-1.137 


-1.419 


-2.200 


-2.134 


CM 

JAN | FEB | MAR | APR , M 




! AUG | SEP 


OCT | NOV 


DEC 






Fig. 29 - 



- Annual equilibrium temperature and thickness field 
for Case 22, maximum snow depth = 120 cm. The scale 
is not suitable for comparison with Fig. 6. 



A 13- 

MAXIMUM SNOW DEPTH = 120 cm. ALL OTHER INPUT DATA AS IN CASE 1. 



J™ Aug Sep Oct Nov Dec Year 



-0.29 


-0.10 


-1.50 


-11.42 


-22.21 


-30.96 


-32.32 


-19.52 


t5.44 


-0.10 


-1.07 


-6.21 


-9.12 


-11.43 


-13.15 


-9.84 


701.9 


702.6 


702.9 


703.0 


703.0 


702.8 


702.3 


702.2 


-0.108 


0 


0.111 


0.306 


0.360 


0.426 


0.394 


2.614 


-0.165 


-0.170 


-0.160 


-0.138 


-0.117 


-0.101 


-0.092 


-1.531 


3.827 


3.750 


2.643 


0.605 


0.085 


0 


0 


16.338 


-1.944 


-0.932 


-0.890 


-0.356 


-0.360 


-0.981 


-1.178 


-15.097 



Date 



Onset 



5 April 
7 June 
12 August 
18 August 
4 October 



Bottom accretion 
Snow ablation 
Ice ablation 



Surface ice ablation 
Bottom accretion 



Other Data 

Surface ice ablation" 0.5 cm 

Bottom ablation 2.1 cm 

Bottom accretion 2.6 cm 
Short-wave radiation penetrating 

surface during snow-free period 0.143 k( 



Equals net bottom accretion. 
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30 -- Average equilibrium thickness of ice as a function of 
maximum annual snow depth. 
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In addition to the previously discussed possibilities, the energy 
balance could be influenced by artificial modification of surface al- 
bedo, long-wave emissivity, incoming long-wave and short-wave radia- 
tion, or the turbulent fluxes. In the remainder of this section we 
shall examine these approaches. 

4.3.1 Influencing the turbulent heat fluxes 

Although the flux of latent heat (F^) and the flux of sensible 

heat (F g ) are small compared with the radiative fluxes, they are not 

negligible when compared with the net energy balance at the surface. 

According to Doronin (see Table 1) , the annual sum of F and F is 

2 si 
only -0.5 kcal/cm year. At first it appears unlikely that this small 

amount of energy could significantly influence ice thickness; however, 
closer scrutiny reveals that variations in F g and F^ can play a large 
part in determining surface conditions, particularly during the abla- 
tion period. Case 23 illustrates these points by arbitrarily assum- 
ing that F s = F £ = 0 during the entire year. The results are given 
in Table 24 and Fig. 31. 

Because the sensible heat flux normally supplies about 0.75 
kcal/cm month to the surface during the winter, surface temperatures 
during this time are 2 to 3°C lower in Case 23 than Case 1; the annual 
average temperature remains unchanged. Over 3 kcal/cm 2 is presently 
lost by the surface during the ablation season by means of the turbu- 
lent fluxes. Eliminating this loss results in 85 cm of ice ablation 
and an equilibrium thickness of 107 cm. This is roughly the same 
thickness change that was achieved when F was increased to 4.5 
kcal/cm year (Case 15) or when a ± is reduced by 10 percent (Case 26). 
The ice would have been even thinner in Case 23 if F and F„ had been 
set to zero during the melt period only. However, in interpreting 
this experiment, it should be borne in mind that the input data for 
the turbulent fluxes in the air boundary layer are not only the least 
reliable ( Fletcher , 1965; Vowinckel and Taylor , 1964), but also the 
most sensitive to those feedback effects that are not accounted for 
in the model . 
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Although there appears to be no direct way of appreciably influ- 
encing F g , the flux of latent heat could be lessened by reducing evap- 
oration during the summer. If evaporation could be completely elimi- 
nated between June and August, an additional 20 to 25 cm of ice abla- 
tion would occur, reducing equilibrium thickness to between 150 and 
175 cm. This suggests that, at least for the Central Arctic, ice 
thickness might more effectively be reduced by retarding evaporation 
from the ice, rather than from the Norwegian Sea. 

4.3.2 Influencing the incoming radiative fluxes 

The current advances being made in weather modification offer 
hope that cloudiness may be controlled to some extent in the future. 
Cloudiness plays the major role in determining the amount of long-wave 
and short-wave radiation reaching the surface. Without specifying 
what cloud conditions would be required, two cases involving increases 
in the incoming radiation fluxes were investigated. 

In Case 24 the incoming short-wave radiation (F ) is made greater 
by 10 percent, or 7.5 kcal/cm 2 year, while all other assumptions are 
unchanged from Case 1. Figure 32 and Table 25 show the results. The 
net short-wave radiation increases by only 1.6 kcal/cm 2 year, because 
of the high albedo of the surface. Snow melt begins earlier and I 
is only 0.4 kcal/cm 2 year greater. Thus the solar energy available 
for melting and heating at the surface is only about 1.2 kcal/cm 2 year 
greater than in Case 1. Ice ablation is 19 cm more at the surface and 
5 cm more at the bottom, resulting in 120 cm less in the average equi- 
librium thickness. The greater surface temperatures are due primarily 
to the lesser ice thickness. 

For Case 25 the incoming long-wave radiation (F^) has been in- 
creased by 10 percent from October through April, simulating a greater 
winter cloudiness than in Case 1. Table 26 and Fig. 33 show the re- 
sults of this integration. Surface temperatures for October through 
April are about 5°C higher, but they remain relatively unchanged during 
the rest of the year. Net long-wave radiation is also the same since 
adjustments in the surface temperatures compensate for a greater F 
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Surface ablation is unchanged, but bottom ablation is greater by 
60 percent. The greater bottom ablation is then responsible for an 
85-cm decrease in equilibrium thickness. 

Although it is not clear to what extent artificial cloud modifi- 
cation could influence F^ or F^, the changes assumed in Cases 24 and 
25 are probably within the range of natural variations. It seems un- 
likely that artificial changes would be much greater. Controlling 
cloudiness thus appears to be one of the less promising avenues to 
ice control. 

4.3.3 Modifying the surface 

Because of the magnitude of the incoming solar radiation, a small 
reduction in surface albedo would provide a large increase in available 
energy at the surface. Dispersal of some dark substance, such as coal 
dust, over the surface might accomplish this reduction. A more feasi- 
ble method might be to introduce snow lichen, which could then propa- 
gate itself naturally ( Fletcher , 1965). In Case 26 it has been assumed 
that, by some unspecified means, the albedo has been reduced by 10 per- 
cent from 1 June until the onset of snow accumulation. Albedos during 
the rest of the year are assumed to be the same as in Case 1. Results 
of the integration are given in Fig. 34 and Table 27. Comparison with 
Case 1 shows an excess of net short-wave radiation of 4 kcal/cm 2 year, 
with an attendant deficiency in average equilibrium thickness of near- 
ly 200 cm. Surface ablation has doubled, bottom ablation tripled, and 
I q is 1 kcal/cm 2 year greater. Temperatures and the dependent fluxes 
behave as would be expected. for ice of 100 cm thickness. 

When summer albedo was decreased by 20 percent (Case 27) , the 
ice vanished during the summer of the third year, before equilibrium 
could be established. 

If the experiments suggested in Cases 26 and 27 were actually 
carried out, it is probable that the resulting thickness changes would 
be even greater than predicted. The dispersal of any dark substance 
on the surface would naturally reduce I q . If I q were reduced by one- 
half, an additional 15 cm of surface ice ablation would occur, enough 
to reduce the equilibrium thickness in Case 26 to less than 50 cm. 
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Table 24. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature ( C) 


-34.01 


-35.00 


-33.25 


-21.97 


-6.08 


Mean ice surface 
temperature ( C) 


-12.43 


-13.53 


-13.79 


-10.70 


-5.05 


Mean ice thickness (cm) 


108.6 


123.0 


136.7 


148.2 


154.8 


Heat flux through surface 
(kcal/cm 2 ) 


1.283 


1.244 


1.094 


0.608 


0.045 


Heat flux through bottom 
(kcal/cm 2 ) 


-1.055 


-1.034 


-0.957 


-0.728 


-0.372 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-1.283 


-1.244 


-1.504 


-2.507 


-3.179 


JAN , FEB , MAR , APR , 




UL | AUG | SEP 


1 OCT | NO 


. . 1. | 






Fig. 31 -- Annual equilibrium temperature and thickness field 



_CASE 23, F g = = 0. ALL OTHER INPUT DATA AS _IN CASE 1. 



-0.10 -0.10 -0.33 -7.79 -19.19 -29.43 -32.01 -18.27 

-0.27 -0.10 -0.28 -4.24 -6.89 -9.34 -11.03" -7.30 

149.0 103.0 66.6 58.4 66.6 79.2 93.9 107.3 

-0.099 -0.187 -0.127 0.693 1.002 1.265 1.280 8.102 

-0.065 0.203 0.340 -0.317 -0.822 -1.026 -1.074 -6.907 

5.678 4.959 2.987 0.605 0.085 0 0 19.740 

-2.015 -0.865 -1.267 -1.297 -1.086 -1.262 -1.281 -18.790 



Date 

1 June 
10 June 
14 June 
26 August 

9 September 



Onset 
Snow ablation 

Ice ablation • 

Bottom accretion 



Bottom accretion 



Surface ice ablation 



Other Data 

Surface ice ablation' 
Bottom ablation 
Bottom accretion 

Short-wave radiation penetrating 
surface during snow-free period 



84.5 cm 
15.0 cm 
99.5 cm 



Equals net bottom accretion. 
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Table 25. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 24, 10% INCREASE IN 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature ( C) 


-30.11 


-32.24 


-30.85 


-21.06 


-7.95 


Mean ice surface 
temperature ( C) 


-13.99 


-15.05 


-15.20 


-12.04 


-6.63 


Mean ice thickness (cm) 


166.7 


177.0 


187.3 


196.7 


202.9 


Heat flux through surface 
(kcal/cm 2 ) 


0.959 


0.998 


0.880 


0.486 


0.060 


Heat flux through bottom 
(kcal/cm 2 ) 


-0.779 


-0.786 


-0.778 


-0.644 


-0.402 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.447 


2.082 


3.442 


Net long-wave radiation 
(kcal/cm 2 ) 


-2.145 


-1.736 


-2.018 


-2.771 


-2.596 


JAN | FEB | MAR | APR , M 






| OCT | NOV 


1 DEC 





100 - 1 




Fig. 32 



-- Annual equilibrium temperature and thickness field 
for Case 24, 10 percent increase in the amount of 
incoming short-wave radiation. 



THE AMOUNT OF INCOMING SHORT-WAVE RADIATION. ALL OTHER INPUT DATA AS IN CASE 1. 



-0.10 -0.10 -0.94 -9.45 -20.02 -28.20 -29.48 -17.54 

-0.95 -0.10 -0.70 -5.20 -9.44 -12.34 -13.52 -8.76 

203.2 170.6 144.0 137.6 138.0 145.9 156.3 168.9 

-0.076 -0.073 0.179 0.757 0.853 0.999 0.973 6.994 

-0.180 0.041 0.150 0.109 -0.460 -0.741 -0.802 -5.271 

5.520 5.455 3.039 0.665 0.094 0 0 20.742 

-2.000 -0.903 -1.037 -0.867 -0.861 -1.552 -1.758 -20.243 



Date Onset End 

2 June Snow ablation 

21 June Ice ablation 

23 June Bottom accretion 

21 August Surface ice ablation 

1 October Bottom accretion 



Other Data 

Surface ice ablation' 58.9 cm 

Bottom ablation 10.7 cm 

Bottom accretion 69.6 cm 
Short-wave radiation penetrating 

surface during snow-free period 1.715 kcj 



Equals net bottom accretion. 



Table 26. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 25, 

THROUGH APRIL. ALL OTHER 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature ( C) 


-26.08 


-28.15 


-26.78 


-17.95 


-8.61 


Mean ice surface 
temperature ( C) 


-13.31 


-14.25 


-14.30 


-11.18 


-6.93 


Mean ice thickness (cm) 


198.5 


206.2 


214.2 


221.6 


226.7 


Heat flux through surface 
(kcal /cm 2 ) 


0.766 


0.812 


0.705 


0.366 


0.081 


Heat flux through bottom 
(kcal /cm 2 ) 


-0.615 


-0.631 


-0.633 


-0.539 


-0.371 


Net short-wave radiation 
(kcal /cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal /cm 2 ) 


-1.951 


-1.551 


-1.803 


-2.460 


-2.303 




Fig. 33 - 



•- Annual equilibrium temperature and thickness field 
for Case 25, 10 percent increase in the amount of 
incoming long-wave radiation during October-April. 



10% INCREASE IN THE AMOUNT OF INCOMING LONG-WAVE RADIATION OCTOBER 
INPUT DATA AS IN CASE 1. 



Aug Sep 



-1.58 -0.10 -0.96 -5.62 -8.09 -11.29 -12.87 -8.37 

229.1 208.7 190.2 183.0 181.1 184.1 190.8 202.8 

-0.083 -0.069 0.213 0.678 0.624 0.794 0.767 5.653 

-0.209 -0.032 0.068 0.083 -0.150 -0.460 -0.603 -4.091 

4.431 4.959 2.710 0.605 0.085 0 0 18.217 

-1.960 -0.903 -0.978 -0.728 -0.624 -1.347 -1.552 -18.160 



Date Onset End 

7 June Snow ablation 

28 June Ice ablation 

1 July Bottom accretion 

18 August Surface ice ablation 

13 October Bottom accretion 



Other Data 

* 

Surface ice ablation 40.6 cm 

Bottom ablation 8.2 cm 

Bottom accretion 48.8 cm 

Short-wave radiation penetrating ^ 
surface during snow-free period 1.291 kcal/cm 



Equals net bottom accretion. 
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Another possibility is that long-wave emissivity (e L > and absorp- 
tivity might also be reduced. Since outgoing long-wave radiation is 
considerably greater than incoming, a decrease in e L would reduce the 
net heat loss from the surface. It is not obvious, however, to what 
extent the resulting temperature increase would affect ice thickness. 
In Case 28, all external parameters are the same as in Case 1, with 
the exception that e L has been specified as 0.8, a reduction of 20 
percent. The equilibrium temperature field is shown in Fig. 35. The 
average surface temperature (Table 28) is 2°C higher than in Case 1, 
and monthly temperatures at the ice/snow interface are up to 7°C high- 
er. Surface ice ablation has increased by one-third, and the average 
equilibrium thickness decreased to slightly more than 50 percent of 
Case 1. The net short-wave radiation has increased by 1 kcal/cm 2 year 
because of the shorter time required to melt the snow and the result- 
ing lower albedo. The net long-wave radiation actually becomes slight- 
ly more negative because of the higher surface temperatures. 

The thickness decrease in Case 28 can be attributed to two fac- 
tors: higher temperatures and increased surface ablation. The in- 
creased ablation is due to greater net short-wave radiation and a 
less negative net long-wave balance during the melt season; the higher 
temperatures are a result of the thinner ice and the smaller amount of 
upward long-wave radiation. Although it is not possible to separate 
these effects, a general idea of the efficiency of each may be obtain- 
ed by a comparison with other cases. The ice ablation and surface 
temperatures in Case 7 (I q = 0) were roughly the same as in Case 28, 
but the thickness decreased by only 45 cm as compared to 130 cm in 
the present case. During the first year that e L was decreased, month- 
ly surface temperatures increased by about 1.5°C in the fall and win- 
ter. Since the ice thickness was still about 275 cm, this increase 
must result from the change in alone and not from a thickness de- 
crease. Thus raising the surface temperature by only 1.5°C caused 
roughly an additional 85 cm decrease in the equilibrium thickness. 
Although we have shown that ice thickness is quite sensitive to sur- 
face temperature, this temperature should be distinguished from the 
air (screen) temperature, about which no statement may be made at this 
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The cases discussed in this section suggest that modification of 
the snow or ice surface is the most effective means of large-scale 
ice removal. Ideally, the surface should be covered with a substance 
that reduces not only the albedo and I , but evaporation and long-wave 
emissivity as well. However, in addition to the logistical problems 
involved in such a'project, it may be difficult to find a material 
with a long-wave emissivity substantially less than one. Furthermore, 
a finely distributed dark solid, like coal dust, would rapidly melt 
into the ice and lose its effectiveness. The ideal material would be 
dark, nontoxic, lighter than water, slowly soluble in water, and have 
a low emissivity. A systematic search should be made to find a sub- 
stance with an optimal combination of these properties. 



Table 27. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 26, 10% REDUCTION 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature ( C) 


-29.41 


-31.58 


-30.40 


-21.32 


-8.77 


Mean ice surface 
temperature ( C) 


-10.86 


-12.08 


-12.45 


-10.01 


-5.62 


Mean ice thickness (cm) 


107.1 


119.5 


131.8 


142.3 


148.8 


Heat flux through surface 
(kcal/cm 2 ) 


1.104 


1.131 


1.010 


0.613 


0.154 


(kcal/cm 2 ) § 


-0.914 


-0.923 


-0.888 


-0.687 


-0.392 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.893 


3.129 


Net long-wave radiation 
(kcal/cm 2 ) 


-2.289 


-1.870 


-2.108 


-2.709 


-2.377 


JAN | FIB | MAR | APR , MA 




1 AUG | SEP 


OCT | NOV 


DEC 






Fig. 34 



-- Annual equilibrium temperature and thickness field 
for Case 26, 10 percent reduction in the surface 
albedo during June-August. 



IN THE SURFACE ALBEDO JUNE THROUGH AUGUST. ALL OTHER INPUT DATA AS IN CASE 1. 



-0.26 -0.10 -0.73 -5.10 -7.07 -8.95 -10.10 -6.94 

141.8 95.7 63.6 59.7 69.1 81.1 94.4 104.6 

-0.103 -0.172 0.098 0.785 0.995 1.163 1.135 7.913 

-0.058 0.240 0.381 -0.524 -0.827 -0.955 -0.964 -6.510 

7.138 5.841 3.194 0.605 0.085 0 0 22.290 

-1.999 -0.903 -1.022 -0.834 -0.994 -1.716 -1.920 -20.742 



Date Onset End 

1 June Snow ablation 

12 June Ice ablation 

12 June Bottom accretion 

21 August Surface ice ablation 

3 September Bottom accretion 



Surface ice ablation 78.3 cm 

Bottom ablation 15.1 cm 

Bottom accretion 93.4 cm 

Short-wave radiation penetrating 2 

surface during snow-free period 2.245 kcal/cm 



Equals net bottom accretion. 
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Table 28. COMPUTED VALUES FOR EQUILIBRIUM CONDITIONS IN CASE 28, 20% 



Variable 


Jan 


Feb 


Mar 


Apr 


May 


Mean snow surface 
temperature ( C) 


-27.97 


-30.31 


-28.93 


-19.34 


-6.98 


Mean ice surface 
temperature ( C) 


-12.78 


-13.85 


-14.01 


-10.93 


-5.92 


Mean ice thickness (cm) 


158.9 


168.6 


178.6 


187.5 


193.2 


Heat flux through surface 
(kcal/cm 2 ) 


0.903 


0.955 


0.839 


0.452 


0.048 


Heat flux through bottom 
(kcal /cm 2 ) 


-0.746 


-0.757 


-0.753 


-0.612 


-0.372 


Net short-wave radiation 
(kcal/cm 2 ) 


0 


0 


0.406 


1.892 


3.129 


Net long -wave radiation 
(kcal/cm 2 ) 


-2.089 


-1.693 


-1.936 


-2.548 


-2.270 


,., JAN i ™ .i MAR „.i, ,. APR \3, 


V 1 JUN | JUL 


AUG | SEP 




DEC 






Fig. 35 -- Annual equilibrium temperature and thickness field 
for Case 28, 20 percent reduction in the long-wave 
emissivity . 



_ 



REDUCTION IN THE LONG-WAVE EMISSIVITY. ALL OTHER INPUT DATA AS IN CASE 1. 



-0.08 -0.10 -0.84 -9.03 -19.25 -26.67 -27.63 -16.43 

-0.67 -0.10 -0.63 -4.96 -9.00 -11.50 -12.47 -8.07 

192.6 161.1 136.5 130.1 130.9 138.9 148.9 160.5 

-0.072 -0.086 0.141 0.734 0.827 0.955 0.924 6.620 

-0.151 0.065 0.164 0.105 -0.480 -0.729 -0.777 -5.045 

5.166 4.959 2.789 0.605 0.085 0 0 19.031 

-1.603 -0.731 -0.841 -0.784 -0.827 -1.509 -1.709 -18.540 



Date 



1 June 
20 June 
20 June 
22 August 
29 September 



Snow ablation 
Ice ablation 



Bottom accretion 



Bottom accretion 
Surface ice ablation 



Other Data 

Surface ice ablation 55.4 cm 

Bottom ablation 11.4 cm 

Bottom accretion 66.8 cm 

Short-wave radiation penetrating 

surface during snow-free period 1.605 kcal/cm 

* 

Equals net bottom accretion. 
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Table 29. SUMMARY OF CASES COMPUTED 
[The "Standard Case" is based on Fletcher's data (Table 1, p. 27). 
For all subsequent cases, only the deviations from 
the input data of Case 1 are listed.] 









Mean Annual Ice 










Thickness 




Case 


Table 


Page 


(cm) 


Characteristics 


1 


2 


54 


288 


Fletcher's heat budget and 

F w - 1-5 kcal/cm 2 year (The "standard case" 


2 






310 


Uniform, low ice salinity (0.09 °/oo) 






66 




Uniform, low ice salinity (0.09 °/oo) 
and F w = 4.5 kcal/cm 2 year 


4 


6 


68 


349 


Uniform, low ice salinity (0.09 %o) and 
near-fresh water (-0.1°C) below the ice 


5 






no ice 


Vowinckel /Orvig heat budget (Table 7 p 71) 


6 


8 


72 


-560 


Vowinckel/Orvig heat budget (Table 7) 
and albedo of Marshunova (Table 1) 


7 


9 


76 


243 


I Q = 07. of F r . No fraction of the short-wave 
radiation penetrating the ice surface 


8 


10 


78 


262 


! 0 = 8.5% of F r during snow-free period 


9 


11 


80 


324 


I Q = 25.5% of F r during snow- free period 


10 


12 


82 


368 


I q = 34.0% of F r during snow-free period 


11 


13 


86 


229 


I o = 34.0% of F r during snow- free period 
and ice albedo lowered to 0,58 


12 


14 


90 


561 


F w = 0. No heat flux in the ocean 




15 


92 


391 


F w - 0.75 kcal/cm 2 year 


14 


16 


94 


162 


F w = 3.0 kcal/cm 2 year 


15 


17 


96 


93 


F w = 4.5 kcal/cm year 


16 






no ice 


F w = 6.0 kcal/cm 2 year 


17 


18 


102 


305 


No snow (h_ = 0). Ice albedo - 0 75 if 
below 272. 9°K surface temperature 










18 


19 


104 


319 


h max = 20 cm 


19 


20 


106 


283 


h = 60 cm 








317 


max 

h max = 80 cm 


21 


22 


110 


411 


h max = 100 cm 


22 


23 


112 


702 


h max = 120 cm 


23 


24 


118 


107 


F s = F £ = 0. No flux of sensible or 
latent heat 


24 


25 


120 


169 


F r increased by 10% 


25 


26 


122 


203 


F L (incoming long-wave) increased by 10% 
during October-April 


26 


27 


126 


105 


Albedo reduced by 0.1 during June-August 


27 






no ice 


Albedo reduced by 0.2 during June-August 


28 


28 


128 1 


160 


Long-wave emissivity reduced from 1.0 to 0.8 
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V. DISCUSSION AND SUGGESTIONS 
FOR FUTURE RESEARCH 

To assist the reader in making his own comparisons among the 28 
cases, Table 29 displays, in briefest form, the parameters varied for 
each case and the resulting equilibrium thickness for each. Also, in- 
side the back cover is an overlay that permits direct visual comparison 
of Case 1 with any other. 

Numerical integration of the sea-ice model has shown that it is 
internally consistent, and, within its inherent limitations, capable 
of predicting the behavior of sea ice under a variety of postulated 
conditions. The experiments described in the previous chapter have 
afforded an insight into the complex relationship between ice thick- 
ness and both the dependent and independent energy fluxes. Besides 
the information generated internally by varying the model, there are 
areas of consistency or inconsistency between the model and field ob- 
servation that merit discussion. The performance of the model has 
been examined for varying values of the input parameters, and the re- 
sults compared with present observations. In this way, it has been 
possible to suggest probable values for many of the external param- 
eters. Analysis of the results has also pointed out features of the 
energy budget that deserve priority in future field studies. 

Wittmann and Schule (1966) believe that open leads are so exten- 
sive in the Arctic that most of the heat transfer between the ocean 
and the atmosphere would occur in these areas. The results of Cases 
1 and 12 to 16 indicate that 1 to 2 kcal/cm 2 year must be supplied to 
the ice by the ocean. It is estimated from sparse oceanographic data 
that only 2.0 to 2.5 kcal/cm 2 year are advected into the central Arc- 
tic Ocean (Crary, 1960). If the model is correct, it is not possible 
to reconcile Wittmann and Schule' s observations and the oceanographic 
information with the present ice thickness. Further study is needed 
to clarify this problem. 

Information about the penetration of the ice by short-wave radia- 
tion (Untersteiner, 1961; Chernigovskii . 1963) is generally deficient. 
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This is unfortunate, as Cases 1 and 7 to 11 have shown the extent to 
which I q can influence ice temperatures, ice thicknesses, and the an- 
nual pattern of mass changes at the boundaries. For a given oceanic 
heat flux, the equilibrium thickness is largely determined by the 
amount of surface ablation, which in turn is strongly influenced by 
the values of I q and the albedo. Although I q was assumed to be con- 
stant during the snow-free period, it is probably a function of time, 
depending upon the changing character of the surface. This time de- 
pendence probably does not affect the annual thickness. 

To understand the interaction of solar energy with the ice, it 
is not sufficient simply to measure incoming and net short-wave radia- 
tion. These measurements should be combined with observations of sur- 
face ice ablation and calculations of I q , preferably by both the in- 
direct method of Untersteiner and the direct method of Chernigovskii. 
Given the uncertainties in the current estimates of average ice thick- 
ness and ablation at the surface, it has not been possible to estab- 
lish definitely the magnitude of I q from the model. On the basis of 
observed and calculated temperature profiles (Figs. 5 and 6) and the 
calculations of Untersteiner (1961), we are inclined to accept a value 
for I q of 2.0 to 2.5 kcal/cm 2 year and a slightly lower ice albedo 
than suggested by Marshunova (Table 1) . 

Representative surface albedos for the ice during the summer ab- 
lation season pose a major problem in treating the arctic energy bud- 
get. We have seen in Cases 5 and 27 that an albedo of 0.45 results 
in a rapid disappearance of the ice. However there is accumulating 
evidence that large-area averages for July may indeed be this low. 
Measurements from high towers (Langleben, 1968) and calculations using 
satellite data (Raschke, unpublished) indicate similar values. How 
is it possible to resolve predictions of the model with these albedos? 
Unless these estimates are all grossly in error, the explanation must 
lie in the way in which melt ponds and leads absorb energy and inter- 
act with the surrounding ice. 

Typical albedos over melting ice range from 0.55 to 0.70; an in- 
tegrated average over bare ice, melt ponds, and leads must be lower. 
In Section 2.4.5 a tentative argument was presented to show that melt 
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ponds have little net effect on equilibrium thickness. If this argu- 
ment is valid, then for a given oceanic heat flux, ice thickness is 
determined by ablation of the bare ice, and hence only ice albedos are 
important in determining equilibrium thickness . Even if the average 
ice thickness is not directly dependent upon the amount of solar en- 
ergy absorbed by the leads and melt ponds, their presence will influ- 
ence the magnitude of the other energy fluxes over the bare ice and 
thus will affect ice ablation to some extent. Presumably this influ- 
ence is included in the albedo values given by Marshunova and Doronin 
(Table 1). 

An important process needing clarification is the way melt ponds 
dispose of absorbed short-wave energy. A melt pond with an albedo of 
0.20 absorbs 17 kcal/cm during July and August, about 10 kcal/cm 2 
more than the bare ice. Since melt ponds are usually about 40 cm in 
depth ( Untersteiner , 1961) , the ponds have to lose roughly 7 kcal/cm 2 
more energy during this period than the bare ice. This energy must 
be lost through larger turbulent fluxes, more outgoing long-wave ra- 
diation, and perhaps most important, a greater I under the melt ponds, 
but the magnitudes of the required increases are unknown. 

Although we have suggested one possible way to interpret the dis- 
crepancies in summer albedos and have pointed out some basic questions 
regarding the melt ponds-, final analysis can only be made on the basis 
of observational evidence. A comprehensive program to define the 
characteristics and behavior of melt ponds is needed. Energy budgets 
over large melt ponds should be determined and compared with similar 
ones for bare ice. The distribution, the depth, and the duration of 
melt ponds should be investigated, along with measurements of temper- 
ature distribution in the water and I q under the ponds. It would also 
be of considerable interest to determine how rapidly the ponds freeze 
in the autumn. If the differences between observed temperatures and 
those predicted by the model are due to a latent-heat release by the 
melt ponds, then it is expected that complete freezing does not occur 
until early December. 

Neither the energy budget of Vowinckel and Orvig (1966, 1967) 
nor that of Fletcher (1965) seems entirely satisfactory. The use of 
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Fletcher's heat budget in Case 1 produced remarkable agreement with 
observations. Paradoxically, it is for this reason that it cannot be 
considered truly representative of a large-scale average. Because of 
the effects of leads and melt ponds, a large-scale average cannot be 
applied to a one-dimensional model of the ice. By assuming an average 
summer surface albedo of 0.64, Fletcher underestimates the average 
amount of short-wave radiation absorbed at the surface in the Central 
Arctic. Although the extent of leads is uncertain, it is known that 
melt ponds cover 20 to 40 percent of the ice during July. Thus, 1 to 
2 kcal/cm 2 year more energy is absorbed at the surface than Fletcher 
assumes. Year-to-year variations appear to be considerably larger 
than this. How this affects his suggested values for the other energy 
fluxes has not been determined. 

The heat budget of Vowinckel and Orvig describes the gross fea- 
tures, but its quantitative accuracy remains to be verified because 
some of the energy fluxes are implicit in their formulation of the 
energy balance. For this reason, their budget cannot be used to de- 
scribe physical changes at the surface. 

The most serious limitation in the present model lies in the sim- 
plifications made in stating the boundary conditions. The independent 
energy fluxes all depend, to some extent, on physical conditions at 
the boundaries. In the present model, only the absorbed short-wave 
radiation and the outgoing long-wave radiation respond to changes at 
the upper surface. The turbulent fluxes, at least, should show a sim- 
ilar dependence. To solve this problem it would be necessary to re- 
formulate the boundary conditions and express the turbulent fluxes in 
terms of the gradients of humidity, temperature, and wind speed in 
the atmospheric boundary layer. If temperature, humidity, and wind 
speed were specified at the top of the boundary layer, it might be 
possible to link the ice model to a model of the atmosphere. A one- 
dimensional model of the ice could then supply realistic information 
at each surface grid point of a model of the atmosphere. Before this 
is attempted, it may be more expedient, as a next step, to develop a 
suitable one-dimensional model of the oceanic convection under the ice 
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and under a free water surface. Combined with the sea-ice model, 
this would enable us to approach the central problem of the transi- 
tion from an ice-free to an ice-covered state of the Arctic Ocean. 
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Appendix 
THE COMPUTER PROGRAM 

The program described in this section is written in FORTRAN IV 
and has been used with the IBM 7090/709 A System at the University of 
Washington Research Computer Laboratory. The program has been gener- 
alized so that the physical parameters of the ice, grid spacing, snow 
accumulation, and the independent energy fluxes are treated as input 
data. Several options allow considerable latitude in the specifica- 
tion of input and output data. The flexible design of the program 
permits it to be easily extended. With only minor modifications, one 
could, for example, include a time-dependent oceanic heat flux or a 
depth-dependent extinction coefficient, or incorporate more complex 
boundary conditions . 

The program is divided into five basic sections: the main pro- 
gram, and Subroutines YARIT, FLIP, SALPR, and RITE. The main program 
reads in the input data, initializes the internal variables of the 
program, and summarizes the results of each year's integration. Sub- 
routine YARIT calculates the temperature and thickness changes of the 
ice and snow for each time step during the year. Subroutine FLIP 
takes the monthly values of the independent energy fluxes at the up- 
per boundary ( F r > F L > F g » and produces smoothed values for each 
time step, according to the method described in Section 3.3. Subrou- 
tine SALPR calculates the salinity profile at each time step [see 
Eq. (34)]. Finally, Subroutine RITE writes the temperature profile, 
ice thickness , and mass changes for each 10-day period throughout the 
year. 

The FORTRAN IV symbols for the input parameters are defined 
below, the units being modified for convenience of machine computation: 

SPH Specific heat of pure ice (cal °K/gm) 

RC Density x specific heat for pure ice (cal/cm^ °K) 
3 

DESS Snow density (gm/cm ) , T < 273°K 
DESJ Density of melting snow (gm/cm^) 
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QS1 Latent heat of fusion for fresh water (cal/cm 3 ) 

QI Latent heat of fusion for pure ice (cal/cm 3 ) 

QB Latent heat of fusion for sea ice at the lower boundary 
(cal/cm 3 ) 

CONO Thermal conductivity for pure ice (cal/cm °K hr) 
EXC Bulk extinction coefficient for sea ice (cm -1 ) 
EXLW Long-wave emissivity and absorptivity of the surface 
WST Melting temperature of snow (°K) 

WIT Melting temperature of sea ice at the upper boundary (°K) 
WTAB Temperature of sea ice at the ice/water interface (°K) 
ABM Average albedo of melting sea ice 

SL1 Ice salinity at the upper boundary [(gm/cm 3 ) • 10 3 ] 

Ice salinity at zj (see Section 3.4) [(gm/cm 3 ) • 10 3 ] 

Ice salinity at z^ (see Section 3.4) [(gm/cm 3 ) • 10 3 ] 

Ice salinity at the lower boundary [(gm/cm 3 ) • 10 3 ] 

DT Time increment (hrs) 

DZ Depth increment (cm) 

Percent of net short-wave radiation which penetrates the 
ice during the snow-free period 



SL2 
SL3 
SL4 



PF 



TMSN Amount of snow which must melt before the snow layer 
is assumed to be isothermal (cm) ? 

GAM Constant ( Y ) in Eq. (3) (cal °K gm) 

BTA Constant (6) in Eq. (4) (cal cm 2 /gm hr) 

SBC Stefan-Boltzmann constant [(cal/°K^ i 

MZZ1 Number of cases to be integrated 
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Snow accumulation (as described in Section 2.4.4) is specified by the 
following symbols : 

SCI Date of Spring increase in the rate of snow accumulation 
(day; e.g. , 1 May = 120) 

SC2 Date of the onset of fall snow accumulation (day) 

SC3 Date of winter decrease in the rate of snow accumula- 
tion (day) 

SC4 Date of the end of snow accumulation (day) 
TSI1 Amount of winter (SC3 to SCI) snow accumulation (cm) 
TSI2 Amount of spring (SCI to SC4) snow accumulation (cm) 
TSI3 Amount of autumn (SC2 to SC3) snow accumulation (cm) 
The independent energy fluxes may be specified by either of two 
methods. If MFX * 0, Subroutine FLIP is called. The fluxes of latent 
heat (Fl) , sensible heat (F2) , incoming long-wave radiation (FDL) , and 
incoming short-wave radiation (FIN), expressed in cal/cm 2 month, are 
read into the program. Smoothed values for FIN and Fl + F2 + FDL are 
then generated for each time increment (DT) and transferred back to 
the main program. If punched-card output is also desired, IPN is set 
equal to zero, otherwise IPN must be set equal to 1. The punched card 
output of FLIP may be read directly into the main program during fu- 
ture runs by specifying MFX =0. In the main program, FSN corresponds 

to F and FLX to F T + F + F„ . 
r L s I 

Termination of program execution is determined by the choice of 
MNY, MXJ, and EQP . The program may be stopped at any predetermined 
time by specifying the year and iteration number. If, for example, 
it is desired that the integration end on 1 July of the tenth year, 
then MNY = 10 and MXJ = number of time increments between 1 January 
and 1 July. Normally, MXJ is set equal to zero. Equilibrium is de- 
fined through EQP. The ice thickness on each successive 1 January 
is compared; when the year-to-year change is less than EQP (0.1 cm 
in the cases discussed in Chapter IV), equilibrium is assumed to have 
been achieved. Before the execution is terminated, the ice thickness 
and temperature profile, along with the internal variables of the pro- 
gram, are stored on logical tape 4. 
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Specif ication of LRT completes the initialization of the program. 
If LRT = 1, logical tape 4 is read. It is thus possible to continue 
the integration of an interrupted case. If LRT = 0, corresponding to 
a new case, it is also necessary to specify: 

ALBD Monthly albedo (ABM overrides the specified values 
during the snow-free period) 

FS Flux of sensible heat (cal/cm^ month) 
2 

FL Flux of latent heat (cal/cm month) 
ITK Initial ice thickness (cm) 

SSTM Initial ice temperature at the snow surface (°K) 

SIT Initial temperature at the ice/snow interface (°K) 

If R5 = 1, temperature profiles will be printed on every tenth 
day. If R5 = 0, only Ice and snow thicknesses and the mass changes 
will be printed. A summary of results is printed out each year, re- 
gardless of the value of R5 . When the ice has attained equilibrium, 
the annual temperature— and-thickness field is printed out for 5 day 
intervals and deposited on logical tape 3 to facilitate numerical 
plotting. 
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MAIN COMPUTER PROGRAM 
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SEA ICE STUDIES-MAIN PROGRAM 
DIMENSION BSAI11, 121 

2ALBD(12),FS(12»,FL(12),T3(600) ' SN,12,60, ' IMM,5 ' 4 >' 

2/BLS2/SLl,SL2.SL3,SL4,Al A2 , P ,u ^ 'uf U rIr, 1 ^Ilr I K t ' VTK ' GEP ' AD0 ' R5 
3ITC,II8,IDT,IZl,IT0,llK,IT4 JJ in 5^ ITT1 * ITT2 ' ITT3 ' IT ^, 

FORMAT ( 12F5 . 0 I 

FCRMATI 12 F6 .0 ) 

FORMAT(12F5.2> 

FORM ATI I4.2F6.1 I 

FC RM AT { 12F6.1J 
) F0RMAT(12(2X,F8.6» » 
> FORMAT! 80H 

' ll^T^n' 70 0C » N - F "-"»« BALANCE FOR ENTI RE ,CE S 

3VERAGE ALBEDO, 2 3X, 8F 8. 2/16H NET sSt" Lvp M^, ?EMP * l 3F8 ' 2/ 1 5H A 
^NETRATING SW , 6X , 12 F8. 1 . F9 3/V 11 m , n f • 5 X , 12 F8 . 1 ,F 9. 3 /15H PE 
5LUX THRU SURFACE »3X,12F8»1 F9 3 /1^7h pP,^ V E»7X f 12F8.0 ,F9. 3/X8H F 
63/21H SURFACE FLUX BALANCE ! 1X*13FS ?/ l^L ArT BOT TOM , 4X , 12F 8. 1 , F9. 
7 R 8.2/19H 3 AVERAGE ICE TH I c" ^ a 3F 8.1 / /, * CTUAL B0T GR0WT H , ,X, 1 3F 

REWIND A 
MZZ = 0 
; READ(5,15) 

CALL FLIP(EXLW) 

GO TO 49 

DO 46 1=1,12 

READI5.6) IFLXI I, j),j=i f60> 
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DO 47 1=1,12 

READ(5,6I ( FSN ( I , J ) , J=l , 60 ) 
DO 50 1=1,12 
DO 50 J=l,60 

FSN< I,J»=FSNII,J>*0T/12.0 
FLX( I,JI=FLX(I,J)*0T/12.0 
SBC=SBC*DT*6XLW 
ITT1=SC1*24.0/DT 
ITT2=SC2*24.0/DT 
ITT3=SC3*24.0/DT 
ITT4=SC4*24.0/DT 
SL1=SL1*0.001 
SL2=SL2*0.001 
SL3=SL3*0.0C1 
SL4=SL4*0.001 
K10=12.0/CT 
Al=<SL4+SLl)/2.0 
A2=l SL4-SL1 >/2.0 
PI=3.1415«265 



Ul=AL0G(0.5«-< (ARSIN( I 2. 0* SL2- SL4-SL1 ) / ( SL4-SL1 ) ) ) I? I ) ) 
U2 = AL0G(0.5 + ( ( ARS I N( ( 2.0*SL3-SL4-SL1 )/< SL4-SL1 ) ) ) /PI | ) 



HAB=0.0 
ShlO=0.0 
L1=0 
ICZ=DZ 
ICT=DT 

TS0P=TSI1*TSI2+TSI3 

ADD1=TSDP*10.0 

ACC=AODi 

ISI=(ADD1/DZ}+1.0 

JJ5=8640/I0T+1 

J17=600/ICZ-1 

J18=ISI+1 

FZ1=10.0/OZ 

IZ1=FZ1 

FI8=360.0/DT 

I I8=F 18 

KK8=II8/3 

FK8=KK8 

ITC=II8*2 

ITO=II8*2/3 

IT4=II8*8 

FW=FW*DT/(24. 0*360.0) 

DSl = TSU*OT/< (SC1 + 360.0-SC3)*24.0I 

0S2=TSI2*DT/C < SC4-SC 1 ) *24. 0 ) 

DS3=TSI3*0T/( (SC3-SC2)*24.0) 

DV1=3.0*SBC 

DV2=C.04*S8C 

XX0=SBC*( (WST/100.0»**4) 
XX1=SBC*< (WIT/100. 01**41 



DES=CESS 

CONS=0.00 68*(DES**2I*3600.0*DT 

CONO=CONO»DT 

BTA=BTA*DT 

A=1.0/(DZ*0Z) 

YY0=FV«/J2.0*QB) 

IF(LRT.EQ.O) GO TO 105 

RfADU) TF,BSM,T1,SL,SND1,RIT,DTB,ABL,EPM,ALBD,IMM,FS,FL 

READC4) NX1,MX2,ITR,L0X,KK1,KK2,JJ1,JJ2,JJ3,JJ4, ILB.ISB.ISF.LCI.K 
15,K7,MMl,fM2,MM3,MIF,Ml,K9,Kll 

READU) 0ZT,DZX,DZP,DZY,EPTC,EPC,DES,RCSP,QS,SOG,FI,FIC,FIOC,TSD, 
1XX2,VA,VB,VC,VD,AIST,ASST,GEP,FTK,TAC «.n»^ 

LCI=1 

Ml = l 

IF{ I7R.LT.JJ5-1) GO TO 100 

ITR=1 

K7 = 0 

GO TO 165 
JA4=ITR 
GO TO 170 

REA0(5,4) (ALBDU ) ,1=1 , 12) 

READ(5,3) JFSU ),I=1,12) 

READ(5,2» (FL(I), 1=1,12) 

REAC(5,5) ITK.SIT.SSTM 

VD=0.0 

AIST=0.0 

ASST=0.0 

DES=DESS 

RCSP=DES*SPH 

ILB= 1 SI+C ITK/IDZ) 

I SB= I S I 

F TK= I TK 

JJ1=ISB<-1 

KK2=ISB-1 

JJ2=ILB-1 

JJ4=JJ2+1 

Ml = l 

OZP=CZ 

DZY=0.0 

SN0=TSI3+ (0S1*C (3 60.0-SC3)*24. 0)/DT) 

I=SND/DZ 

KK1=ISB-I-»1 

ISF=KKl-2 

ZZ = I 

OZX=SND-ZZ*DZ 
DZ7=DZX+DZ 
XX2=0ZX 
MX1=10 

CALL SALPF(FTK,MXl,TACJ 
Y=(WTAB-S IT )/FTK 
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X=0.0 

DO 110 I=ISB, ILB 

Til I l=SIT*X*Y 
110 X=X+OZ 

CO 135 I=!LB,600 
135 T1(I)=WTAB 

DO 140 1=2, ISF 
140 T1(I-1)=H£T 

TKISF) = SSTM/100.0 

C0N=C0NO+ ( { BTA* SL { JJ 1 ) l/ITl! JJ1 1-273. 0> » 
Ti(KK2l=Tl( ISB)*(CCN*(T1(ISB)-T1UJH) > /CONS 
VA=VA+ALBC(LOX) 
AP=CCNS/( SND-OZ) 

BP=- ( CONS*T 1( KK2 ) I /( SND-DZ )-FLX< 1 ,1 > 
145 X=T1(ISF) 

TH I SF 1 = 1 CV1*(X**4 )-BP)/((DV2*(X**3 H-AP)*100.0) 
Y=ABS(X-TH ISF) » 
IFCY.GT..C000C01J GO TO 145 
V0=VD+FLXtl,l)-(SBC*ITH ISF)**4I > 
T1USF)=T1USF)*1CC.0 

Tl( ISF+1) = T 1 C ISF)*((DZX*CT1 (KK2 l-Tl ( I SF I ) )/(SNO-DZ) ) 

A1ST = AIST*T1( ISB» 

ASST = ASST4T1( ISF) 

Y=(T1(KK2)-T1( ISF-t-1) )/{ SNO-DZT) 

X=0.0 

DO 150 I=ISF,KK2 

T1(I*1)=T1( ISF>1)>X*Y 
150 X=X+DZ 

CCN=CONO*C18TA*SL(JJ4) » / t T 1 (J J4 )-273. 0 ) » 

Tl(JJ2l=WTAB-( (FW*DZ)/CON) 
152 TF<1,1)=0.0 

DO 155 K=2,60 

I = ( (K-l )*10»/IDZ+i 
155 TF( 1,K )=213.0— Tl( I ) 

DC 160 K=l,12 
160 ABL(K)=O.C 

K5 = l 

K7 = l 

EPTC=C.O 

ITR=2 

FIC=0.0 

GEP=0.0 

FIOC=0.0 

MIF=0 

MN1 = 0 

MM 2=0 

MM3=1 

VA=0.0 

VB=0.0 

VC=0.0 
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161 QS=QS1*DES 

SOG=CONS/ (DZ*DZ*RCSP> 

T SD=SND 

SNDH1» = TJD 

RIT< 1» = FTMGEP 

DTB ( 1 »=FTK*ADD+GEP 

ABL ( 1 ) = 0. C 

EPM( 1»=0.C 

LCI=1 
165 MX2=4 

MX1=10 

R1K=RIT(1 » 

L0X=1 

JM=1 

FI=0.0 

EPC=0.0 

K9=0 

Kll = 0 

VTK=0.0 

170 CALL YARIT(MX1,TSI3,TMSN) 

TSDP = TSIHTSI2*TSI3 

ADD=ADD1 

J 20= I SB- I SI 

I SB= ISI 

KK2=ISI-1 

JJ 1= I SI+1 

ISF= ISF-J20 

KK1=KK1-J20 

ILB= ILB-J20 

JJ2=ILB-2 

JJ3=ILB-3 

JJ4= ILB-1 

DO 190 1= ISF, ILB 

J=I+J20 
190 T1(I)=T1(J) 

DO 1S5 !=ILB,J17 
195 Tl(It=WTAE 

DO 200 I=I5B,ILB 

J=I*J20 
200 SL(I)=SL(J> 

*RITEt6,10> (SLm,I=ISB,ILB) 

WRITE(6,18) LCI 

hRITE(6,15) 

CO 205 1=1,5 

X= IMM C I t 1 ) 

J=X/«2.0*FI8» 

IMHI I,2)=J+1 

Y=J 

Z={X-2.0* Y*FI8)*DT/24.0 
IMM( I,3)=Z 
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Y= IMM I » 3 ) 
205 IMM(I,4)=CZ-Y)*24.0 

IFCLQX.EQ.13) GO TG 212 

DO 210 K=L0X,13 

00 210 J=l,ll 
210 BSM(J,K)=0.0 

GG TO 220 

212 DO 214 1*1,11 
X=0.0 

DO 213 J=l,12 

213 X=X+BSM(I,J) 

214 BSMI,13)=X 
BSM1,13)=BSM(1,13)/12.0 
BSMJ 2,13I=BSM(2, 13 1/12.0 
BSK 11,13 )=BSM( 11, 131/12.0 
00 215 1=2,8 

215 BSM(I,13)=BSM(I,13)/100C.O 

22C WRITE (6, 17) LC I , I I ( 1 , I ) , I =2 , 4 ) , ( I MM( 2 , 1 ) , 1 = 2 , 4 ) , t I Mf { 3, I ) , 1 = 2 , 4 1 
1, ( IMM( 4 , I ), 1 = 2,4) , ( I MH{ 5, I )»I=2,4) 

WRITE (6, IS) «BSM( 1,1 ) ,1 = 1,13) , (BSMI2, I ) ,1=1 ,13) , (BSMI3, I ), 1=3,10), 
1 (BSM (4, I ), I =1,13), (BSM(5,I 1,1=1,13) ,(BSM(6, I 1,1=1,13) ,(8SM( 7,1), 1 = 
21, 13),(BSM(8,I ),I=1, 13) , IBSM19.I) ,1=1,13) ,(BSMJ 10,1 ) , 1= 1 , 1 3 ) , { BSM( 
311,1 ),I = i,13) 

X=BSM(5,13)-BSM(8, 13 )-BSM ( 7 , 13 ) - ( SWLO/1000. 0 ) 

WRITE(6,16) SWLO.X 

SWLO=0.0 

WRITE(6,15) 

LCI=LCI*1 

IFIMNY.LE .LCI ) GO TO 225 

ITR=1 

K7=0 

IF(ABS(RTK-RIT< 1 I > .GE.EQP) GO TO 165 

MXJ = 0 

R5=1.0 

MNY=LCI 

TF(1,1)=0.0 

DO 227 1=2,60 

J=( ( 1-1 )*10)/IDZ*l 
227 TFC1, 11 = 273. O-TKJ) 

GO TO 165 
225 IF(MXJ.EQ.O) GO TC 230 

MXJ = 0 

ITR=1 

K7=0 

R5=l.0 

GO TO 165 
230 CCNTINUE 
235 REWIND 4 

WRITE (4) TF ,BSM, Tl ,SL, SND1 ,R I T , DT 8 , A8L , EPM , ALBD , I MM , F S, FL 
WRITE14) fXl,MX2,ITR,L0X,KKl ,KK2,JJ1 ,JJ2, JJ3.JJ4, ILB , ISB , I SF , LC I , < 
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15,K7,MM1, MM2,*M3,MF,M1 ,K9,K11 

WRITE (4) CZT f DZX,DZP,DZY,EPTC,EPC,DES,RCSP,QS,SOG,FI , FI C ,F I OC.TSD. 
lXX2,VA f VB,VC,V0,AIST,ASST,GEP,FTK,TAC 

MZZ=MZZ+1 

IF(MZZ.LT.MZZl) GO TO 925 
REWIND 3 
STOP 
END 



YEARLY ITERATION SUBROUTINE 
SUBROUTINE YAR I T ( * XI , T S I 3 , TMSN ) 
DIMENSION T(600» , TT(600I 

COMMON /BLKl/TFl 12,80 ),BSM (11,13 ) , Tl ( 6001 , EPM 1 1 2 J , A8L ( 1 2 ) , S ND1 ( 12 » 

I. DTB(12>,RIT(12>/BLSl/SL(6C0l/BLAl/FLXC12,60),FSNU2,6OI ,IMM(5,4), 
2ALB0(12),FS(12),FL(12),T31600> 

CEMMON /BLSK/ISB,IDZ,DZ,DZX/BLK2/ISF,ISI,KK1,JJ2,JJ4, J17, J18,M1 ,LC 

II, ITP.MXI »K5,DZP,0ZY,EPC,EPTC,TSD,TSDP,WTAB,WST,FTK,VTK,GEP,ADD,R5 
2/BLS2/SLl,SL2,SL3,SL4,Al,A2,PI,Ul,U2,ILB/BLA2/ITTl,ITT2,ITT3,ITT4, 
3ITC, 1 18, I CT, 121 ,IT0,ITK,IT4, JJ1 , JJ3 , JA4 , J.J5 ,K7 , KK2 , KK8, LOX , LRT , MNY 
4,^XJ,MX2,fMl,MM2,MM3,MIF,NR,ABM,AIST,ASST,A,Ll,;;TW , A BI , BT A ,CONO, C 
50NS,SWLO,DT,DESS,DESJ,DES,DVl,DV2,DSl,DS2,DS3,DZT, DTI, EXCEPT, EPS, 
6FI,FIC,FIC,FW,FI0C,FI8,FK8,FZ1,GAM,PF,QI,QB,QS,QS1,RC,RCSP,RTK, SPH 
7, SOG.SBC, SN0,TAC,VA,VB,VC,VD,WIT,XX1,XX5,XX2,XX0, YYO , K9 , K10 ,K11 

FORMAT! 16 F8.3 > 

F0RMAT(3I5,10X,2F10.2) 

FCRMAK15H SHEARING ERROR ) 

00 565 JJ=JAA,JJ5 

IF< ITC.NE.K7) GO TO 100 

BSM(1,L0X » = 273.0-(ASST/(2.0*FI8) ) 

BSM(2,L0X > = 273.0-(AIST/(2.0*FI8» I 

BSM(3,LOX)=VA/(2.0*FI8) 

BSMU.LOX ) = VB + FICC 

BSM(5,L0X )=FIGC 

B SMI 6, LOX )=VD-FS( LOX)-FL(LCX) 

BSM( 7, LOX )=VC 

BSHl 8, LOX l = FIC 

BSM<9,LOX)=BSM(4,L0XI-BSM<5,LOXH-BSM(6,LOXJ*BSMC7,LOX)+FS(LOX)+FL( 
1L0XI+VMF 
VMF=0.0 

BSM{ 10,LOXI=EPC-FI 

BSM( ll,LOX)=VTK/(2.0*FI8) 

VTK=0.0 

VA=0.0 

VB=C.O 

VC=0.0 

VC=O.C 

FIOC=0.0 

ASST=0.0 

AIST=C.O 

FIC=C.O 

LCX=L0X+1 

K9 = 0 

K11 = 0 

IF(12.LT.L0XJ GO TO 570 

FI=EPC 

K7 = 0 

K7=K7+1 

KS=K9+1 

IFU9.LE.K10.AND.K9.NE.1 ) GO TO 104 
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K9 = l 

K11=K11+1 
FLW=FLX<LOX f Kll ) 
F£W=FSN(L0X,K11 ) 

10 IF( ITT4.LT.ITR) GC TO 115 

TSD=TSD+DS2 

XX2=XX2+DS2 

GO TO 180 
■5 DZX=0.0 

DZT=DZ 
MX1=1 
MX2=1 
TSD=TSCP 
KK1=ISF+1 
GO TO 200 
0 IF( ITT2.GE. ITR) GO TO 185 
MX1=6 
DZX=0.0 

ALBD(L0XI=ALB0(L0X+1| 

MX2 = 2 

MIF=0 

TSD=0.0 

JJI=KK1+1 

CES=OESS 

QS=QS1*CES 

PCSP=DES*SPH 

C0NS=0.OOea*(DES**2)*3600.O*DT 
SOG=CONS/ <OZ*DZ*RCSP) 
> IF( I TT3.L T.ITR ) GO TO 150 
TSD=TS0+0S3 
CZX=DZX+DS3 

IFJ2.0*DZ.LT.TS0) GO TO 140 
IEIDZ.LT.TSDI GO TO 13-= 
ISF=IS8-1 
GO TO 380 

IFIMX2.EQ.3) GO TO 380 

ISF=ISB-2 

MX2 = 3 

DZX=DZX-D* 
GO TO 380 
MX2 = 4 
MX1 = 8 
ISF=IS6-3 
XX2=TSD-2.0*DZ 
KK1= IS8-1 
KK2=ISB-1 
JJ1=IS8+1 
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145 DZX=XX2 

0ZT=DZ+XX2 

GO TO 200 
150 MX1=7 
155 TS0=TSD*DS1 

DZX=CZX*DS1 

GO TO 130 
160 IF( ITT3.LT. ITR> GC TO 165 

TS0=TS0+DS3 

XX2=XX2*DS3 

GO TO 180 
165 MX1=9 
170 1SD=TS0+0S1 
175 XX2=XX2*DS1 
180 IF(CZ.GE.XX2) GO TO 145 

XX2=XX2-0Z 

KK1=KK1-1 

ISF=ISF-1 

GO TO 180 
185 X=( 1.0-A8M>*FSW 

VC=VD+FLW 

FIO=X*PF 

FI0C=FIOC*FI0 

FNW=FLW-F IO + X 

VA=VA+ABM 

VB = VB + X-F 10 

190 RCP=RC+lG/M*SUKKl>)/( ( T 1 ( KK1 I -273. 0 » **2 > 
C0N=CQNO+ ( ( BT A* SL ( KK1 ) )/(Tl(KKl J-273.0) ) 
SIG=A*CON/RCP 
AP=CON/( ( 1.0*SIG)*0Z) 

BP=-FNW-< CON*( ( 1.0-SIG)*T1(KK1I+SIG*T1 (KK1+1) ) ) / ( ( 1 . 0 +S I G ) *DZ ) 
Y=T1I ISI ) /ICO.O 
195 ET=(DV1*( Y**4)-BP)/( ( D V2* ( Y** 3 » + AP ) * 1 00. 0 ) 
X = AB S ( ET— Y ) 
Y = ET 

IF(X.GT..C0CC0C1) GO TO 195 
VO=VD-( S8(* (Y**4 ) ) 
T( ISI ) = Y* 100.0 
IFIWIT.GE .T( ISB ) ) GO TO 380 
T( ISI> = WIT 
GO TO 380 
200 X=I1.0-ALeD(L0X) )*FSW 
VC=VD+FLW 
VA=VA+AL8D(L0X) 
V6=VB+X 
FNW=FLW*X 

IFtCZX.EQ.O.OJ GO TO 205 
AP=CONS/( (1.0+S0G>*DZ+DZX| 

BP=-FNW-( CCNS*( ( 1.0-S0G»*T1(KK1 )+SOG*Tl (K.KH-1 ) ) ) / ( ( 1 . C+ SOG ) *DZ*-DZX 
1 I 
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GO TO 210 
205 AP=CCNS/( (1.0*S0G)*DZ) 

BP=-FNW-( CONS*( <1.0-SOG(*T1(KK1 ) «-SOG*Tl ( KKH-1 ) ) )/ ( ( 1 . C+SOG ) *DZ ) 
210 Y=T1( ISF ) /100.0 

215 ET=( DV1*( Y**4)-BP) /( (DV2*( Y**3 )+AP> *1 00. 0 ) 
X=ABS(ET-V) 

Y = E T 

IFIX.GT..C00C001) GO TO 215 
T( ISF)=Y*10C.O 
IF(WST.GE.T( ISF) ) GO TO 219 
VC=VD-XXO 
M>1=1 
MX2=1 
GO TO 235 
219 GG TO(225,38C,38C,220) ,MX2 

2 2C TCKKl-l >=( ( 1.0*SGG»*DZ*T (I SF)+DZX*( ( 1 .0— SOG )*T1 IKK1 )+SOG*Tl(KKl*l) 
1 ) »/( (1.0 + SOG)*OZ*DZX) 
T3(KK1-1I=T(KK1-1) 

225 00 23C K = KK 1 1 KK 2 

230 T(K)=U l.C-SOGt*TKKH-SOG*(TUK + l)+T(K-l) ))/(1.0+S0G) 

VC=VD-(SBC*((T( I SF )/100.0)**4) ) 

GC TC 380 
235 IF( MM1.NE .0) GO TC 240 

MM1 = 1 

IMM(1,1)=1TR 

IFITSD.GE.TSCP) GC TO 240 

TS0P=TSO 
240 T(ISFI=WS1 

GC T0(270,245,255) ,MX1 
245 IF(KK1.GT.KK2) GO TO 260 

CO 250 K = KK 1 » KK2 

T(K)=WST 
250 TT(K)=WST 
255 T(ISB)=WI7 

TTt ISB)=WIT 

X=ABM+(TSC*ABI/TSCP» 

VC=VD+FLW-XXO 

Y = < 1.0-XMFSW 
FNW=FLW*Y 
VA=VA+X 
VB=VB+Y 

EPT = (XX0-FNV»)/QS 
VMF=VMF+EPT*QS 
GC TO 275 
260 MX1 = 3 

GO TO 255 

265 RCP = RC4(GAM*SL(KK1 »)/« (TKKKl )-273.0)**2) 
CON = CONO*(< BTA* SL ( KK1 1 > / ( T 1 ( KK1 ) -2 73. 01 ) 
SIG=A*CON/RCP 
T( ISB)=WIT 



-153- 



X=( 1.0-ABF»*FSW 

VC=VD+FLW-XXl 

F IO = X*PF 

FIGC = FIOC-»FIO 

FNW=FLW-FIO+X 

VA=VA+ABM 

VB=VB+X-F 10 

X = (XX1-FNU/(2.0*QI) 

Y = (( 1.0*S1G>*CZ+DZX)/2.0 

EPT=X-Y+SQRT( U + Y)**2 + C0N *(WIT*( SI G-l. 0 ) *T 1 ( KK1 ) - S IG*T1 ( KK 1+1 > ) /Q 
II I 

VPF=VMF+EFT*QI 
GO TO 275 
270 X=(XXC-FNH/(2.0*CS) 

Y=I I 1.0 + SCG )*0Z*DZX)/2.0 

EPT=X-Y+SCRT( (X+YI**2+CCNS*(WST+< SOG-1.0)*T1 ( KK 1 ) -SOG "»T 1 ( KKU 1 1 )/U 
IS) 

VMF=VMF+EFT*QS 
275 EPTC = EPTOEPT 
TSD=TSD+EPT 
CZT=CZT+EPT 
CZX = CZT-DZ 

IF1MX1.GE.2I GO TC 285 
IF(EPTC.GT.TMSN) GO TO 285 
MX1=2 
DES=OESJ 

TSD=CESS*1SD/DESJ 

ABI={ <ALBC(LOX)-AEM)*TSCP»/TSD 

EPTC=TSD-TSDP 

K=TS0/DZ 

Z = K 

KK1= ISB-K+1 

ISF=ISB 

DZX=TSD-CZ*Z 

0ZT=CZ*DZ> 

QS=QS1*DES 

RCSP = DES*SPH 

CCNS=C.OOfc8*( DES** 2 1*3600. 0*DT 
SCG=CONS/ (DZ*DZ*RCSP ( 
HI«1 = 0 

DO 280 K=KK1,KK2 
T { K I =ViST 
280 TT(K»=WST 
T(ISB)=WIT 
T T ( I SB I =W IT 

285 IFIOZT. GE.OZ.AND.CZT.LT. 2. 0*0Z) GO TO 305 
IF(DZT.GE.2.0*DZ) GO TO 300 
CZT=CZT+02 
KK1=KK1+1 
CZX=DZT-DZ 
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GO TO(295 ,265, 285,2901 ,MX1 
290 JJ1=KK1+1 

ISB=KKl-2 

Ll=l 

ISF=ISB 

GO TG 285 
295 IF(KK1-ISF-2.LE.C> GO TO 285 

ISF=ISF+1 

GO TO 285 
300 DZ1=C2T-D2 

KK1=KK1-1 

ISF=ISF-1 

CZX=CZT-DZ 

GC TO 285 
305 GO T0C315 , 375, 310,330) ,MX1 
310 IFITSD.GT.O. ) GO TO 375 

MX1=4 

I SB=KKl-2 

ISF=ISB 

M IF=1 

IMM( 2,1I=ITR 
JJ1=KK1+1 
GC TO 330 

315 T(KK1 ) = ( ( 1.0-SOG»*T11KK1)+SOG*(T1(KK1*1 )♦ (1.0-DZX/DZT )*WST) 1/(1.0+ 
1S0G*<1.0-CZX/DZT» » 

T7(KK1)=T(KK1I 

IFIEPTC.GE.O. > GO TO 320 

T(KKl-l)=VST+( ( (T(KK1)-WST»*DZX)/DZT) 

T3(KKl-l >=T(KK1-1) 
320 J20=KK1+1 

DO 325 K=J20,KK2 
32 5 T(K)=( ( 1.C-S0G>*T1(K)*S0G*(T1(K*1)+TIK-1> ) )/U.0+SOG> 

GC TO 375 

330 RCP=RC*( G*M*St(KKl )) /( < T I ( KK1 1-273. 0 ) **2 ) 
C0N=CON0 ♦ ( ( BTA*SL ( KK1 ) ) / ( Tl ( KK1 (-273.01 ) 
S1G=A*C0N/RCP 

T(KK1) = < { 1.0-SIG»*U<KK1) + SIG*T1(KK1 + 1)+SIG*(1.0-DZX/CZT)*WIT>/<1. 
10+SIG*l l.C-DZX/DZT) > 
TT (KK 1 1 = T (KK1 ) 

T IKK 1-1 ) = UT+(DZX* (T(KKII-WIT) I /DZT 

T3(KK1-1)=T(KK1-1> 

T3( ISB)=WIT 

I F ( E PT.L T .0. ) GQ TO 375 

MX1=5 

Ll = l 

IMMC 3,1 )= ITR 

VD=VD+XX1 

EPTC=EPTC-EPT 

VPF=VMF-EPT*QI 

ABL( K5 + 1 )=TSDP*EPTC 
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CZT = CZT-EFT 

OZX=CZT-DZ 

K=JJ2 

J2C=KKl-2 

Y=CZ-DZX 

IF<DZX.GT.5.0) GO TO 335 
ISF=KK1-1 
ISB=ISF 
GC TO 340 
5 KKI=KK1-1 
0 T1(ISF)=W!T 

Z=FTK4GEP+TS0P+EPTC 
J=Z/OZ 

ZZ = (J-D* IDZ 

OZP=Z-ZZ 

JJ2=J+ISF-1 

ACO=( ISF-I)*ICZ 

JJ4=JJ2*1 

RB=JJ2 + 2 

JJ3=JJ2-1 

JJI=KK1 

DZX=C.O 

OZT=CZ 

CZY=OZP-OZ 

TAC=DZY 

CO 345 J=KK1,JJ2 
J20=J20*1 • 

I!^"T l ^ 20 '* ,Y * ,T1,J20 * l »- T1 «J201M/DZ 
TllJJ4+l)=w TAB 
Tl( ISB-1 )=hST 
IF(J2C+1-K| 355,35C,360 

^Q^S'lIs tJ2 °* 1,+ ( Y * tT1(J20+ 2»-Tl(J20 + l» , »/DZ 

WRITE<6,12) 

GO TC 365 

IF(J20.NE.K) GO TC 355 

Jo J 3?c =T jiKK^J 1 l2 {Y * ,WTAB " T1,J20 + 1,n/(DZP - DZ ' 
T1(JI=T(J( 

FTK=FTK+GEP+TSDP*EPTC 

CALL SALPPIFTK.MXl.TAC) 

GEP=0.0 

EPTC^C.O 

TSD=0.0 

IF( ITT2.GE. ITR) GC TO 190 

XMTT3-ITR 

DS3=TSI3/X 

VA=VA-ABM 

Ve = VB-»FIO«( (ABM-1.0)*FSW( 
VC=VD-FLW 
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FIQC=FIOC-FIO 

GG TO 120 
75 CCNTINLE 
30 JJ3=JJ2-1 

rnM = ^ MG ' M * SLUJ2n/nT HJJ2)-273.0)**2I 

Sci!icr sLijj2,,/,Ti, - , - i '' 

Y=(( 1.0+SlG)*DZ+DZY)/2.0 
XYT=(T1(J J2 )*WTA8)/2.0 
XYS=CSUJJ2) + SL4)/2.0 
CCN=CONO+(BTA*XYS/(XYT-273.0) I 

CZP=CZP+EFS 
DZY=DZP-DZ 

iG^ j2 D zYr* ,(i -°- si ^ 

IF( l.NE.M IF | GO TC 383 
Z=FTK+GEP-DZP 
CZP=CZP-EPS 
IF(MX1.EQ.5) GO TC 381 
Z=Z+TSOP+£PTC 
1 v C tt^U? AM * SL,JJ2,/,(TT(JJ2, - 27 3.0(**2I) 

xyt--:xS;b^: 0 fio * exp, - 2+exc,/rcp » 

CCN=CCNC+(BTA*XYS/(XYT-273.0» ) 

czp s : C z7^ s ^^ 

CZY=DZP-DZ 

FIC=FIC+(CON*(X-hTAB)/DZP» 
GC TO 384 

^£ = FKMCGN*(TT<JJ2I-WTAB)/DZP) 
GEP=GEP+EPS 

T T Tu^^T T ;L + i! ozYMTTuj2 '- wTABn/ ° zp ' 

IF(EPSI38J,«10,4CC 
IFCNM2.NE.0J GO TC 390 
MM2=1 
M*3 = 0 

IKM(4 f 1 )=ITR 

IF(DZP.GE.DZ) GC TO 410 

CZP=CZP+D2 

CZY=CZP-DZ 

T(JJ2)=TT(JJ2) 

JJ2=JJ2-1 

JJ3=JJ2-1 

JJ4=JJ4-1 

TT(JJ2I=(CZP*T(JJ2 + 1,-DZ*WTAB)/DZY 
IF( ILB-JJ4-1.LE.C > GO TO 355 
ILB= ILB-1 
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5 T( IL B ) = In TAB 

TT ( I L B ) =H TAB 

GO TO 385 
3 IF(MM2.NE.O) GO TC 405 

M*2=0 

MM3=1 

lt*m 5,1 )= 1TR 
5 IF(DZP.LT.2.0*DZ> GO TO 410 
OZP=CZP-D2 
DZ¥=CZP-DZ 
JJ2=JJ2*1 
JJ4=JJ4+1 
ILB=ILB+1 

T( JJ4 )=Tt JJ2 | ♦{ (DZ*UTAB-T(JJ2) H/DZPI 

TTUJ4)=T(JJ4) 

GO TO 400 

3 IFIMX1.GT.1. AND.MX1.LT. 6) GO TO 450 
DO 415 J=JJ1,JJ3 
K=JJ3+JJ1-J 

RCA=RC + (GAM*Sl.(K) (/( ( T 1 ( K ) -2 73. 0 ) **2 ) 
CAN=CONO+((BTA*SL(KI ) / ( T 1 ( K 1-273. 01 ) 
SAG= A*CAN /RCA 

' TT(K)=( ( 1.0-SAG)*TKK) + SAG*lTltK-l)*TTlK+l» I ) / ( 1.0+SAG) 
IF(MX2.EQ.1.0R.HX2.EQ.4I GO TO 435 
RCP = RC*(.GAM*SL(KKl)l/((Tl(KKl)-273.0)**2) 
C0N=C0N0+UBTA*SL(KK1) ) / ( T 1 1 KK1 )-2 73 . 01 ) 
SIG=A*C0N/RCP 
MX3=1 

xis:;^ 

AP=CGN/(DZ*(1.0*XX5) I 

X=( l.O-ALBD(LOX) )*FSW 

VC = VO + F|_W 

VA=VA+ALBC(LOX> 

VB=VB+X 

FNW=FLm-X 

BP=-FNH-C *P*TT(KK1 ) ) 

Y = T1< ISFJ/100.0 

ET = (0Vi*O**4)-BP)/( (DV2*< Y**3>+AP)*100.0) 
X=ABSCET-Y) 

Y = ET 

IFIX.GT..C000001) GO TO 420 
TTIISF)=Y*100.0 
IF{MX3.EQ.2) GO TG 430 

T( ISB»=(TT( ISF>+XX5*TT(KK1))/U.0«-XX5) 
DO 425 J=KK1,JJ3 

RCP=RC+(GAM*SU J I >/((Tl( J 1-273.01 **2) 
CCN=CONO+UBTA*St(J) ) / ( T 1 1 J ) -2 73. 0 H 
SIG=A*CON/RCP 

T«J>»((1.C-SIG»*T1IJ|*SIG»(T1IJ + 1)*HJ-1)|)/C1.0*SIG» 
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T3(KKl)=(1(KKl)+TT(KKl))/2.0 
CON=CCNO«-((BTA*SLCKK1> ) / ( T 1 (KK1 1-273.01 I 
AP=C0N/(0Z*(1.0+XX5) ) 
BP=-FNW-(*P*T3(KK1 M 
Y=TT { ISF ) 
MX3=2 
GO TO 420 
430 T(ISF»=TT(ISF> 

VC=V0-( SBOM! T< ISF)/100.0)**4) ) 
RCP=RC-MGAM*SL(KK1 H/( ( T 1 ( KK1 >-273. 0 ) **2 ) 
CON=CONO+UBTA*SL(KK1H/(T1CKK1 )-273.0) ) 
SIG=A*CCN/RCP 

T3(ISBI = M(ISF|+XX5*T3{KK1))/<1.0+XX5» 
IFIMX2.LT.3) GO TO 460 

T3IISB-1J =T(ISF)+( If T3I ISBI-TC I SFI )*(TSD-DZ) J/TSD) 
GO TC 460 

435 CCN = CONO+<(eTA*SL(JJin/<TllJJl 1-273.0) ) 

T( ISB)=( ( C0N*TT(JJ1) l/CCNS+T(KK2) )/( I . 0+C ON/CON S I 
TT( ISBI=Ti ISB) 
OC 44C J=JJ1,JJ3 

RCP=RC+CG/M*SU J I )/< [Til J J-273. 0 > **2 I 
C0N=CONO+((BTA*SL ( J) )/(Tl( J)-273.0) ) 
SIG=A*CCN/RCP 

44C TU> = (( 1.0-SIGI*T1IJ»«-SIG*IT1(J + 1I*TIJ-1M »/(1.0*SIGI 
CC 445 K=KKl,KK2 
J=KK1+KK2-K 

445 » a «tl.O-SOG)*Tl( J »*S0G*CT1(J-1) + TT(J*1> I)/ II. O+SCG) 

GO TO 460 
450 OC 455 J=JJ1,JJ3 

K=JJ2+JJ1-J 

RCA=RC + ( GAM*SL( K) »/{ f T 1 ( K ) -2 73. 0 > **2 > 
RCP = RC+(GAM*SU J) )/< (TKJ)-273. 0)**2) 
CCN=C0N0*<( BTA*SLU> l/ITH J 1-273. Q) ) 
C AN=CONO + < ( BTA* SI ( M I / ( Tl ( K )-273. 0) ) 

SIG=A*CON/RCP 
SAG= A*C AN /RC A 

TT(K» = ((1^-SAG)*T1{K)>SAG*(T1(K-1)*TT(K + 1)))/ ( 1.0^SAG) 
*55 T(J» = ((1.0-SIG»*T1(JUSIG*(THJ + 1I+T(J-1( ) ) / ( 1 . 0* S I G ) 
460 RCP=RC*fG*M*SL(JJ2*l)»/||TlIJJ3 + l)-273.o!**2 

CCN=CONO*I(BTA*SL(JJ3*1M/(T1<JJ3*1)-273.0> ) 

SIG=A*C0N/RCP " ~ ' ' 

^n J ^s 1 ! = ^J* ( !" 1 , SIG, * T1UJ3+1KSIG * (T1,JJ3+2,+T,JJ3),,/(1 - 0+ SIG) 

LU 46D J=KK1,JJ4 
465 T3(J)=(T( J)+TT< J) )/2.0 

IF(MX1.GT.1.AN0.MX1.LT.8) GO TO 470 

CCN = CCNO+((ETA*SL( JJ 1 1 1/lTllJjl t-273.0) I 
,-,„ I 3 ! ISB,= t <C0N*T3UJ1) )/C0NS*T3(KK2l)/(1.0 + C0N/C0NS» 
470 IFI1.NE.MIF) GO TO 490 

IFIMX1.NE.5) GO TC 475 

Z = Cl 
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J20=KK1 

6C TO 480 
5 Z=DZX 

J20=KK1-1 
0 DO 465 J»J20,JJ4 

PCP-RC'MGAM*SHJ) |/( (T3< J)-273.0)**2) 

T3U >=T3{ J)+{EXC*FIO*EXP«-Z*EXC) ) /RCP 
5 Z=Z*CZ 

ShLO=SWLO*F IC*EXP(-Z*EXC> 
0 DO 495 K=KKl,JJ4 
5 T1(K.| = T3(K) 

GO T0( 500, 520, 520, 51 5, 520,500, 500, 51 5, 51 5, 515, 5 15 I. MX1 
0 GO T0f505,515,51C, 515»,MX2 
5 IF(fcST.GT.TUSFJ) GO TO 520 

T1(KK1-1)*T(KK1-1) 

GO TO 520 
0 Tl(KKl-2)=T3(KKl-2) 
5 TKKK1-1)*T3<KK1-1) 

:• IF(HX1.GT.1.AN0.MXI.LT.5» GO TO 525 

THISF| = T(ISF) 

GC TO 530 
> TUISBI=WIT 
) TAC=TAC+EFS 

IF(MX1.NE.4) GO TC 535 

TAC=TAC-»EFT 

MB=TSDP*EPTC 

GO TO 540 

IF(TAC.LT.OZ.ANO.TAG.GT.O.O) GO TO 545 
HTK=FTK+GEP+HAB 
L1 = 0 
l-AB = 0.0 

CALL SALPP(HTK,MXl,TACI 
ASST = ASST«T1< ISFt 
IFIR5.NE. 1.0) GO TO 555 
IF(MOD( ITR.KK8) > 55 5, 550, 555 
WRITE(2I T1,ISF,ILB,0ZX,DZY 
WRITE(6,1C) IT1UI ,J=1,ILB) 
WRITE(6,U) ISF,ILB,ITR,OZX,DZY 
AIST=AIST«Tll I SB I 
VTK=VTK*F1K-»GEP 

GO TO ( 560 ,558 » 559 , 560 ) , MX2 

VC=VC + ( (CCN*ITHKK1J-THISF)»)/(DZ*DZXI » 

GO TO 563 

VC=VC+(C0NS*(T1(ISBI-T1(ISF))/DZX) 
GC TO 563 

rn = Tn + c^ NS * ,T1<I£B, - T1<ISF »'/<DZ+DZX>| 
GO TO 563 

CCN=CONS 

GO TO 557 
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V7K = VTK*TSDP-f£PTC 
GC TO 557 

CALL RITECMXl.nia, ' 5 " 

GC TC(565,57C),MX3 

ITR=ITR+1 

RETURN 

END 
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C FLUX INTE PPOLAT ICN SUBROUTINE 

SUBROUTINE FLIPIEXLWI 

CCMMON /BLA1/FLXC12,60) ,FSM12,60I , I MM (5,41 , ALBO ( 12 ) , FS ( 12 ) ,FL( 12) 
1,T3<600) 

CIMENSION FIX(60),FSH20) , FLMT ( 60 ) , Fl ( 20 ) , F2 ( 20 ) ,FDL ( 2C ) , FI M 20 ) , F 
H60),FLX2UC),F3(6C) 
2 FORMAT! 12F5.0) 
FORM AT ( F5 .1 , 13 ) 

5 FORMAT! 12F6.1 ) 

6 FORMAT! 10F7.0) 

7 FORMAT! 12 f 6* C ) 

10 FCRMAT! 10JF8.1, 5X1 ) 

11 FORMAT ( 1H 1 , 20X , 9H DELTA T=,F5.1,6H HOURS//) 

12 FORMAT ( 1H C ) 

13 FORM AT ( 1 5 KO TOTAL FL X ( I ) =F 6 . 1 , 10X , 12 H TOTAL F(I) = F8.1) 
READ!5,4t DT.IPN 

REA0(5,7) (Fill 1,1=2,17) 
READ15.2) !F2U >, 1=2, 17) 
READ(5,6) {FDL(I),I=2,17) 
READ(5,6) (FIN( I),I=2,17) 
F 19=720. O/DT 
I I9=FI9 

WPITE(6,11) DT 

LIX = C 

Fill 1=0. 0 

F2( 1 )=0.0 

FCL!1)=0.C 

FIN(1)=0.C 

DO 100 1=2,17 

F1(I)=F1( II/FI9 

F2( I )=F2( II/FI9 

FCL! I )=(FCL( I l*EXLW)/FI9 
100 FIN( I)=FIM D/FI9 

DO 1C5 1=2,17 
105 FKII=FU I)+F2U) + FDL(I( 
X=«!DT/72C.C)**3>/2.0 
X2=720.0/CT 
X3=2.0*72C.0/CT 
X4=3. 0*720. O/DT 
Ml = l 
L1=0 
M4=l 

19=1 19/2 
Kl = 2 
K2 = 5 
ITR = 0 

ue-i 19+1 

DO 2C0 1=2,13 
DC 110 K=l,17 
110 FS1(K)=F1(K) 
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DO 115 K=l,60 
115 FIX(K>=O.C 

120 Al = ( FSH I+l)-FSiU*2) + (FSi(I + 3)-FSlC I»>/3.0>*X 

A2=« ( (X2+>3+X4 1*FSl( I »- I X2+X3 >*F SI ( 1+3 ) ) /3. O* ( X2+X4 I * FS 1 1 1 + 2 I- 1 X3+ 
1X4)*FSH 1*1) )*X 

A3=(X3*X4*FS1( I + n-X2*X4*FSiU+2)+(X2*X3*FSl< I + 3>-(X2*X3+X2*X4+X3* 
1X4>*FS1( I H/3.0)*X 
A4=< X2*X3*X4*FSH I >*X 1/3.0 
Y = FI9 

DO 135 J=1,II9 

L1=L1+1 

Y=Y+1.0 

FLMT (J ) =A 1* ( Y**3 M A2* ( Y**2 ) +A3* Y+ A4 

GO T0(130,125),M1 
125 F ( L 1 ) = FLM K J ) 

GC TO 135 
130 FIX(L1»=FLMT(J) 
135 CONTINUE 

L1=0 

GO T0( 140,150, Ml 
140 CO 145 K=K1,K2 
145 F5HKI = FIMKI 

Ml = 2 

GC TO 120 
150 K1=K1+1 
K2=K2+1 
MR=ITR + 1 
Ml = l 

GG T0U55,U5>,M4 
155 DC 160 J=l, 119 

FLX2CJI=FlX(Jt 
160 F3(J)=F(J) 

M4=2 

GG TC 200 
165 DO 170 J= 1, I 9 
K=I9+J 

FLMTI J)=F3(K) 

FL MT ( K ) =F I J ) 

F3(K l = F1 K ) 

F ( J > = FLMT (J ) 

F ( K I =FL MT (K ) 

FLMTI J»=FLX2(KI 

FLMT (K I =F IX ( J ) 

FLX2(K)=F IX(K) 

FIX( JI = FLPTU» 
170 FIX(K>=FL*T(K) 
175 LIX=LIX+1 

DO 1 80 J= 1, I IS 

FLX(LIX,J ) = FIX( Jl 
180 FSN< LIX.J l = FIJI 
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SUMX=0.0 
SLMF=0.0 
CO 185 J=l, 119 
SLMX=SUMX4FIX( J) 
185 SUMF=SUMF+F(J I 

WRITE (6, 12) SUMX,SLMF 
hRITE(6,12) 

WRITE(6,1C) CFIX(K),K=1,II9) 
WRITE(6,12» 

WRITE<6,iC) »F(K) ,K=1,II9) 
IF(12-ITR)2CC,19C,200 
190 ITR=ITR*1 

DO 195 J=1,I9 
K=I9+J 

F IX ( J ) = FL >2 I K ) 

FIX(K)=FL>2(J) 

F( Jl =F3(K > 
195 FJKI=F3(J) 

GO TO 175 
200 CONTINUE 

CO 2C5 1=1,60 

FSNC 1,1 ) = C.O 

FSN( 2,1 )=C.O 

FSN( 11, I )=0.0 
205 FSN{ 12, I ) =0.0 

IF( IFN.EQ.l) GO TO 220 

DC 210 1=1,12 
210 PUNCH 5, ( FLX( I, J ) , J=l,60) 

DO 215 1=1,12 
215 PUNCH 5, ( FSNI I, J) ,J=1,60) 
220 RETURN 

END 
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C SALINITY PROFILE ADJUSTMENT 

SUBROUTINE SALPR t THK ,MX1 ,TAC,HK1) 
COMMON /BLS1/SL16C0J 

CCMMON /BLSK/ISB,IDZ,DZ,DZX/BLS2/SL1,SL2,SL3,SL4,A1,A2,PI,U1,U2,IL 
1 B 

Zl = ( 10.0*THK)/34.C 
Z2=( 20.0*THK)/34.C 
SL( ISBMSL1 
SL( ILB)=SLA 
EP1=L1/ALCG(Z1/THK) 
EP2=U2/ALCG{Z2/THK ) 

GI = U1.0/EP2-1.0/EP1I*THK»/CZ2-Z1) 

G2=1.0/EP1-(G1*Z1/THK) 

J1=ISB+1 

J2=ILB-1 

IFIMX1.NE.4 ) GO TO 101 

Z = OZX 

GC TO 102 

101 Z=DZ 

102 DO ICO I=J1,J2 
ZP=Z/THK 

EP = 1.0/(G 1*ZP*G2I 

SL ( I )=A1 + /)2*SIN(PI*( (ZP)**EP-0.5) » 

100 Z=Z*DZ 

I=THK/DZ 
Z=I*ICZ 
TAC= ThK- Z 
RETURN 
END 
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C OUTPUT SUEROUTINE 

SUBROUTINE R ITE ( MX 1 , MX3 ) 

COMMON /BLK1/TFU2,80» ,B SMI 11, 13 J ,T1 (600) ,EPM( 12) ,ABL (12) ,SND1( 12) 

I, CTB(12),RIT(12) 

CCMMCN /BISK/ 1 SB, I DZ , DZ , DZ X/BLK2/ I SF.ISI ,KK1 ,JJ2,JJ4, J17, J18.M1.LC 
1 I, ITR.MXl ,K5,DZP,DZY , EPC , E PTC , TSD, TSDP , WTAB ,WST , FTK, VTK , GE P , ADD, R5 

10 FORMATt 13X.2H 1.7X.3H 10.7X.3H 20.8X.2H 1.7X.3H 10.7X.3H 20.8X.2H 

II. 7X.3H 1C.7X.3H 20.8X.2H 1.7X.3H 10.7X.3H 20//I 

11 F0RMAT(2X,I3,5H CP. , 12 ( F8. 3, 2X) ) 

12 FORMATI//10H ACCRE T ION , 1 2 ( F8. 3, 2 X ) / 1 OH BOT OEPTH, 12 ( F8. 3.2X )// ) 

13 FGRMAT{ IH1, 115X.6H YEAR I2/10H SNO DEPTH , 12 ( F8. 3 , 2X I / 10H ICE THICK 
1.12IF8.3, 2X)/10H ICE A8LTN , 12 < F8 .3 , 2X ) // ) 

14 FORMAT! 8H DEP TH , 13 X , 7H JA NUARY ,23X , 8HF EBRUAR Y , 23 X, 5 H PARCH , 25X , 5HA 
1PR IL ) 

15 FCRMAT(8H DEPTH, 15X, 3HPA Y »27X » 4HJUNE » 26 X, 4HJUL Y ,24X ,6 H AUGUST ) 

16 FORM AT ( 8H DEP TH , 12X, 9HSE PTEMB ER , 22 X, 7H0CT0B ER ,23X , 8HN0VEMBER, 22X 
1 , 8HDECEMB ER ) 

17 FORMATl 7H 1 CATE.10X.10H SNO DE PTH , 10 X, 10H ICE THIC K , 10 X, 10H ICE A 
1BLTN.10X, 10H ACCRETION, 10X.10H BOT DEPTH) 

18 F0RMAT(2X ,2I3,9X,5 (F8.3.12X) ) 
K5=K5+1 

EFM( K5 I =EFC 

CTB( K5)=F1K«GEP*ACD 

SNCKK5 ) = 1SD 

IF(MXl.NE.A) GO TC 130 

ABL(K5)=TSDP*EPTC 

SND1(K5) = C0 

RIT(K5I=FTK+A8L(K5 l+GEP 

GO TO 135 
130 RIT(K5)=CIB(K5)-ACD 
135 DO 14C K=JJ4,J17 
140 T1(K+1)=WTAB 

CO 145 K=Z, ISF 
145 T1(K-1)=WST 

IFIR5.EG.C.0) GO TO 155 

DO 150 K=2,6C 

I9=( (K-1)*1C)/IDZ+1 
15C TF(K5,KI=273.0-T1(I9) 
155 I F( I TP. GE .MX I ) GO TO 160 

IF(K5-12)200,16C,160 
160 IF(R5.EQ.C.O) GO TO 191 

WRITE (6, 12) LCI,(SND1(JI,J = 1,12),(RIT(J),J=1,12),UBL(J),J=1,K5) 

GO T0( 170 ,175,180, 1651 , Ml 
165 Ml=l 
170 WRITE(6,1M 

GC TO 185 
175 WRITE(6,15) 

GC TO 185 
180 WRITE(6,1<» 
185 WRITE(6,1C) 
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CC 190 J=l,50 
I6=( J-l)*10 

190 WRITE(6,11> 16, ( TF ( K , J t , K= 1 ,K5 ) 
WRITE(6,12| (EPM(J ) , J = l , 12 I , ( DTB (J ) ,J=1,K5) 
GG TO 194 

191 IFtMl.NE. ] ) GO TO 192 
J6 = l 

J7=-10 
WRITE (6, 17) 

192 CO 193 K=l,12 
J7=J7+10 

IFU7.NE.20) GO TC 193 
J6=J6 + 1 
J7 = C 

193 tcJIf'f: 18 ' J6 ' J7 ' SNDUK »' R IT(K),ABL(K»,EPM(K),DTB(KI 
IFIMl.NE..-) GO TO 194 

M1=0 

194 M1=M1+1 

CC 195 K=],K5 

195 ABL(K)=O.C 

DO 196 1=1,12 
SN01 ( I ) = 0.0 
R I T( I 1=0. C 

196 EPN(I)=O.C 
K5=0 

IF( ITR.GE .HX1 » GO TO 2C5 
200 MX3=1 

GO TO 210 
205 MX3=2 

ITR=ITR+1 
210 RETURN 

END 
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